DEPARTMENT OF OCEAN ENGINEERING 

MASSACHUSEHS INSTITUTE OF TECHNOLOGY 
CAMBRIDGE. MASSACHUSETTS 02139 



AUTOMATIC CONTROL OF A SUBMERSIBLE 



Kurt James Harris 



SM (NAME) 



Course XIIIA 



Thesis 

H2897 



'June 1984 



/ 



f 






AUTOMATIC CONTROL OF A SUBMERSIBLE 



by 

Kurt James Harris 

II 

B.S., Chemistry, Miami University 

Submitted in Partial Fulfillment 
of the Requirements for the 
Degree of 

MASTER OF SCIENCE IN NAVAL ARCHITECTURE 
AND MARINE ENGINEERING 

and the Degree of 

MASTER OF SCIENCE IN MECHANICAL ENGINEERING 

at the 

MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
May 1984 

@ Kurt James Harris, 1984 



AUTOMATIC CONTROL OF A 
SUBMERSIBLE 



by 



Kurt James Harris 



Submitted to the Departments of Ocean Engineering and Mechanical 
Engineering on May 1984 in partial fulfillment of the requirements for 
the Degree of Master of Science. 



ABSTRACT 



Multivariable control theory is used to design a controller for a 
modified drooping stern turning submersible. The procedure used is 
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controller successfully minimizes roll and pitch for a turning submer- 
sible. 
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CHAPTER 1 



INTRODUCTION 



1 . 1 Introduction 

In the present age of high technology and sophisticated develop- 
ment of vehicles, both air and submersible, designers continue to push 
the performance limits of their designs, due to increased demands on 
system requirements* Aircraft are required to perform faster, at higher 
and lower altitudes. Submersibles , especially submarines, are being 
designed to dive deeper and faster. It is this last category of ve- 
hicles that this thesis addresses. 

The operating envelope of a submarine is predicated on its design 
depth and speed plus its ability to handle a mechanical casualty with 
enough time to recover before exceeding test depth of the hull. Cur- 
rently, submarines have the capability to operate at depths greater than 
700 feet and speeds in excess of 30 knots. 

The current hull design has been determined to be the most effi- 
cient hydrodynamic design and has evolved from the Albacore. Addition- 
ally, the rudder and sternplane cnfiguration has basically been un- 
changed from the cruciform design for decades. The cruciform stern is 
quite a sensible design in case of a casualty situation. When one wants 
to climb or dive, all one has to do is either command rise or dive on 
the sternp lanes. To turn, a similar situation exists. One just needs 
to turn the rudder in the desired direction one wants to go which also 
causes the boat to dive. However, the cruciform stern arrangement has 
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its drawbacks. There is no differential control, i.e., both the upper 
and lower rudders are mechanically tied together as are the right and 
left sternplane control surfaces. A severe casualty such as a stern- 
plane jam on dive does not allow for independent action on the stern- 
plane control surface. The corrective action to combat this casualty is 
to quickly reverse propulsion power direction and command full rise on 
the sailplanes. All of this corrective action must be accomplished fast 
enough, before the boat exceeds test depth. Finally, the cruciform 
stern does not effectively allow the boat to maintain a level attitude 
in a turn. The boat has a natural tendency to roll thereby causing the 
rudder to contribute to pitch as if it were a stern surface. 

In order to effectively correct the drawbacks of the cruciform 
stern, an alternative stern configuration could be used employing dif- 
ferential control of its surfaces, i.e., each surface independently 
operated. The independent action would allow the operator to deflect 
the control surfaces that effect depth in a casualty. Also, if a stem- 
plane jam casualty occurred, then the plane that did not suffer the 
casualty could be used to control the rate of depth excursion, while the 
propulsion plant is used to slow the vehicle. Differential action of 
the stern control surfaces can also be used to counteract the forces 
that cause roll, thereby prohibiting undesirable depth changes of the 
vehicle . 

The stern configuration that will be considered in this thesis is 
called a modified drooping stern (MDS). The MDS consists of a single 
rudder surface with two dihedral surfaces 37.5"^ below the horizontal. 

In theory the MDS should be able to effectively control roll and enable 
the operator to combat a sternplane casualty. However, the MDS presents 
the operator with a coordination problem in that its movement is not so 
straightforward; for example, for normal cruising, how should one de- 
flect the dihedral surfaces to control the effects of roll? It there- 
fore becomes of considerable importance to design an automatic control 
system that handles such a vehicle in operation. 
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In this thesis a truly multivariable automatic controller will be 
designed for the MDS configuration. Since there exists a number of 
operating points within the submerged envelope, we have selected as a 
nominal design point that which describes the situation of a two degree 
deflection of the rudder and small deflections of the dihedrals, while 
the boat is at a velocity of 21 knots. 

A sketch of the proposed submersible is located in Figure 1. The 
overall length is 180 feet with a top speed of 21 knots. Mass proper- 
ties are summarized in Appendix A. The convention for positive deflec- 
tions of the control surfaces are shown in Figure 2, using the right 
hand rule system. The rudder and dihedral surfaces will be referred to 
as Dl, D2 , and D3 and will be limited to a maximum deflection of plus or 
minus forty degrees. 

The mathematical model used for the control design utilizes the 
six degree-of-f reedom nonlinear differential equations as developed by 
Gertler and Hagen known as the NSRDC Model 2510 [1]. The linear hydro- 
dynamic coefficients employed in the 2510 equations are derived in Chap- 
ter 2. All nonlinear coefficients used are based on submersibles of 
similar size and shape. Since this submersible has never been con- 
structed on scale basis or tested in a tow tank, the accuracy of the 
coefficients could not be verified. However, the computer model behaved 
in a reasonable manner as one would expect with a full scale equivalent 
model. 

The origin of the submarine fixed axes is located at the center 
of gravity. The positive directions of axes are: x-direction, forward; 

y-direction to starboard; z-direction downward. The linear and angular 
velocity components, forces and moments are shown in Figure 2 [6]. 

This thesis will develop and design a multivariable automatic 
controller for the submarine model in Chapter 2. Chapter 2 consists of 
the development of the linear hydrodynamic coefficients and explains the 
modifications made to the six degree-of-freedcm equations of motion. 
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SCALE: 1/16 in. = 1.875 ft 




11 



t 



M,q 







Figure 2. Sketch showing positive directions of axes, 
angles, velocities, forces, and moments. 
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Additionally, an explanation is provided as to how the equations of 
motion were linearized. Chapter 3 consists of the controller (compen- 
sator) design. It includes a brief explanation of the control varia- 
bles, system augmentation, and calculation of the control gain matrix, 

G, and the filter gain matrix, H, which result from the frequency domain 
loop shaping methodology employed for the design. The design method- 
ology itself is also discussed briefly. Chapter 4 summarizes the re- 
sults of the nonlinear controller simulations. Finally, Chapter 5 lists 
the conclusions and recommendations as drawn from the resultant sub- 
marine model and controller design. 

1. 2 Contributions of Thesis 

The contributions of this thesis is twofold: (i) evaluation of 

the MDS configuration for a submersible in a high speed turn, 

(ii) design of an automatic control system to facilitate the MDS opera- 
tion using. Linear Quadratic Gaussian/Loop Transfer Recovery (LQG/LTR) 
methods . 

The emphasis is on the frequency-domain properties of the state- 
space based LQG designs. Simple and efficient tools for loop shaping 
are employed that lead to satisfactory designs. These tools accomplish 
singular value ’’pinching" in the high and low frequency ranges, so as to 
satisfy stability and performance requirements. 

Moreover, this thesis tests the controller designed on the basis 
of a linear model on the nonlinear process model which is a better 
approximation to the actual system. 
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CHAPTER 2 



THE SUBMARINE MODEL 



2 . 1 Nonlinear Equations of Motion 

The equations of motion describing the design vehicle in this 
thesis as already mentioned, were the six degree-of-f reedom equations 
contained in the NSRDC document 2510 [1]. The coefficients of the equa- 
tions, as developed in Section 2.3, describe the vehicle’s geometry. 
Modification was made to the equations to account for the geometry of 
the MDS configuration. The terms that involved both rudder and stern- 
plane action were combined; these can be seen in the lateral force, 
pitching moment and yawing moment equations. The factors that premul- 
tiply the corresponding surface deflections were based upon the geometry 
depicted in Figure 3. 



D1 is primarily used for rudder control while control surfaces D2 
and D3 are used for depth and roll control. However, D2 and D3 can also 
act as an additional rudder and their potential contribution can be seen 
from the lateral force equation where the multiplying coefficient is 
given by sin 37.5° = 0.6088. The resulting term looks like the follow- 
ing, D1 - 0.6088D2 + 0.6088D3. Additionally 3 this term takes into ac- 
count the sign of the hydrodynamic coefficients and In this 

case both and Y^^ are positive and employ the sign convention as 
established in Figure 2. Likewise, the stern contribution terms in the 
pitching and yawing moment equations have similar form. 
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For depth control, surfaces D2 and D3 are used. The factor 
cos 37.5® = 0.7934 appears. Depending on the sign of the hydrodynamic 
coefficients N^, and the sign of the deflecting 

term is assigned. 



One must keep in mind that^ not all hydrodynamic coefficients 
involving a control surface have a value that can be predicted or meas- 
ured in model tests. If a value could be determined, then one would 
assign the additive or subtractive nature to the contribution of the 
control surface based upon the sign of the corresponding hydrodynamic 
coefficient. The resulting nonlinear and auxiliary force equations are 
given in Appendix B. 



2. 2 Linearized Equations of Motion 

Control design requires that the nonlinear equations of motion be 
linearized about a nominal operating point. The nonlinear equations are 
of the form 



Ex = f (x, ij) 



( 2 . 1 ) 



from which one can derive the linearized equations in the form 

( 2 . 2 ) 

which represents a Taylor Series expansion of the nonlinear equations 
about a nominal point ^ and by retaining only the first order 

terms. More specifically, expanding Equation (2.1) in a first order 
Taylor Series around a nominal point, one gets 



9f (x 



+ EAx = 






u ) + 






9x 



(x - 



8f(x , u ) 

X ) + — - (u - u,) 



9u 



(2.3) 
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If one lets Ax - (x - x ) and Au = (u - u ) where x is the nominal 

--0 — — — o — o 

state vector and the nominal input which, furthermore, are known to 
satisfy the original nonlinear equations, then one is left with 

9f(x , u ) 3f(x , u ) 

EAx = ^ Ax + ^ Au (2.4) 

dx — 8u, — 

where A = r and B = . The linearized equations are 

— d2£ — dU 

listed in Appendix C. 

Additionally, the full nonlinear equations include crossflow and 
vortex shedding terms. For the purpose of this thesis, the crossflow 
term was linearized while the vortex shedding term was omitted. The 
reason for this was to be able to compare the results of this submer- 
sible to a full scale model, whose linearized equations did not include 
vortex shedding, either. 

2 . 3 Hydrodynamic Coefficients 

The hydrodynamic coefficients that comprise the linear and non- 
linear equations of motion are very important in the mathematical model. 
The linear coefficients can be estimated using equations based on the 
geometry of the submersible. The nonlinear coefficients are normally 
the result of model tests either in a tow tank or wind tunnel. The 
coefficients not calculated in this thesis are taken from past designs 
of similar size and shape vehicles. 

2.3.1 Body Coefficients [2] 

We assume a prolate ellipsoid for the hull shape, where the major 
and minor axes are given by 
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a = 90 ft b = 10 ft 

u u2 

V = -j Trab 

V = 37,699 ft^ 



The buoyant force of the ellipsoid is 



B = pV 

B = 2.4283 X 10^ Ibf. 



The surface area of the ellipsoid is given by 



c o . o ab • “1 

S = 2irb + ZTT — sin e 

e 



where 



. . c 

e = eccentricity = 

c 
e 

S = 8932.4 ft^ 



(a^ - b^)2 



0.9938 



From Lamb [3] with ^ “ 9, the factors 
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K' = 0.8645 

= 0.0241 

= 0.9538 



can be obtained. The moments of inertia for the body are as follows 



for a prolate ellipsoid 



m = 



B 



32.17 ^ 
s 



2.4283 X 10 Ibf 

32.17 ^ 
s 



.= 75483 slugs 



2 2 
m(a 4- 



yy 



1.2379 X 10^ slugs-ft^ 



I = 1.2379 X 10® slugs-ft^ 

zz ^ 



m(2b ) 



XX 



I = 3.0193 X 10® slugs-ft^ 

XX ^ 



The mass and moments of inertia can be nondimensionalized by 
dividing with the appropriate factor. 
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Therefore , 



I’ = I’ = 0.000654 

yy zz 



I’ = 0.000016 

XX 



m’ = 0.1293 



The center of gravity is assumed to be located in the port- 
starboard plane of symmetry, 12 inches below the centerline. 



= 0 

Zg = 1 ft 



Nondimensionalizing, 



= ^ = 0.00556 



From Newman [4], the Reynolds number for a body in 10°C seawater at 
21 KTS, is given by 



R = — 



Y 



(21) (1.689) (180) 
1.35 X 10~^ (3.281)^ 



= 4.3931 X 10 
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Assuming that the frictional drag coefficient is equal to the drag coef- 
ficient and using the ITTC frictional resistance line gives 






0.075 



(log^Q R - 2)- 



= 0.0017 

Jj r 



for the coefficient of drag. 

The body hydrodynamic coefficients can now be calculated [5] 



= -0.234(m')°*^^ - = -0.008009 (2.5) 



Y' 

r 


= -(0.10 - K^m' = 0.000981 


(2.6) 


Z' 


= Y' = -0.008009 


(2.7) 


w 


V 




Z' 


= -Y' = 0.000981 


(2.8) 


q 


r 




M' 

w 


= 0.87(K2 - K^m' = 0.010458 


(2.9) 


M' 


= 0.045(m') = 0.000582 


(2.10) 



q 
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-0.010458 



( 2 . 11 ) 



N' = -M' = 

V w 



= -0.000582 (2.12) 

X^’ = -K^m' = -0.000312 (2.13) 

Y^' = -K^m' = -0.012333 (2.14) 

Y- ’Z 

Y-’ = = -0.000069 (2.15) 

= Y-' = -0.000069 (2.16) 

V p 

Z-’ = Y-’ = -0.012333 (2.17) 

w V 

M.» = -K^I’ = -0.000565 (2.18) 

q y 

N-’ = M-’ = -0.000565 (2.19) 

r q 



2.3.2 Fin Coefficients [6] 

The MDS stern configuration introduces asymmetry about the cen- 
terline of the vertical plane. Several views have been depicted to 
allow calculation of areas, aspect ratios and planes of motion. Fig- 
ures 3 and 4 depict fin dimensions. 

For a single fin: 



= Projected area 
= 100.4 ft^ 
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! 



,1 




4.5 ft 





Figure 3. Rear view, rudder, and dihedral fin dimensions. 
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SCALE: 1/8 in. = 1 ft 




PROJECTED AREA = 194 ft^ 

PROJECTED AREA OF LOWER 
PLUME PROJECTION = 28.8 



4.5 ft 



VERTICAL 




Figure 4. Vertical and horizontal plane projections. 
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fi 



A 

m 


Projected area of movable portion 
80.6 ft^ 


AR^ = 


aspect ratio of fin = — 

c 


s = 


28 ft c = 8.25 ft 


AR^ = 


3.394 


AR 

m 


aspect ratio of the movable portion 
6.222 



For the horizontal plane projection contributions, one has: 



= 


projected area of the horizontal plane 
186.4 ft^ 


^ = 


aspect ratio in the horizontal plane 
2.69 



The following definitions are used in calculating A and 



‘"La 


2ir 



In the vertical plane 



ARf = 


3.394 


^f = 


100.4 ft^ 


*"La 


3.953 
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X 

cp 



distance from center of buoyancy to center of pressure 
of the stern 



The center of pressure is assumed to be 1/2 the mean chord length 
forward of the hailing edge. For this vehicle 



X 

cp 



-94.5 ft 



Y' 

V 






-0.012249 



( 2 . 20 ) 



V PP 

N' = = 0.006431 

V i 



( 2 . 21 ) 



The contribution to roll due to asymmetry is 



K’ 

V 



Y’Z 
V cp 

I 



where was calculated by summing the different components of the 

vertical projection times their moment arm and then divided by the total 
area. 



For this vehicle a value for Z was found as 

cp 



Z = 10.9 ft 

cp 



Therefore 



K’ = -0.000742 

V 



( 2 . 22 ) 
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I 

[i 



I 

I 



( 



For the horizontal plane projection 



whereby 



In 

idual fin 



a = 2.64 

A = 186.4 ft^ 

C, = 3.604 

La 



Z' 

w 



M' 



w 



Z’ 

q 



Y' 

r 



K' 



r 



M' 

q 



N' 

r 



-C^ A 
La 



w cp 
Z 



M' = 
w 



N' 



V 



-K' X 
V cp 




££ 



Z 



Y'x 
r cp 



= -0.020733 

= -0.010885 

-0.008673 
0.006276 

= -0.000380 

= -0.004553 

= -0.003295 



order to calculate K' , the method will be applied to an 

P 



(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

indiv- 
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I 

I 



I 

I 



-0.012249 



A = 100.4 



then 



= 3-953 

La 



ytt 

V 



-C^ A 
La 




Z’ is the distance from centerline to center of pressure of the 
cp 

individual fin. For this fin Z’ = 7.44 ft 

cp 



K’ = -0.000021 

P 



The K’ for the design vehicle will be three times the above value. 



K’ = -0.000063 (2.30) 

P 



Abkowitz [7] provides a means for calculating the acceleration 
derivatives. In the horizontal plane 

2 2 

2. = TTpS C 

“ ’ 4^PT7 

where s is the span and c is the chord. 

s = 11.1 ft + 11.1 ft = 22.2 ft c = 8.25 ft 
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t 



0.000331 



(2.31) 



Z- 

w 



In the vertical plane, 



Z- 

w 




I 



Y- ' 

V 



2 2 
TTPS C 




+ C 




where 

s = 14 ft + 8.5 ft = 22.5 ft c = 8.25 ft 

Hence , 

Y*' = -0.000336 

V 



Z* 'x 

W I 



M* ' 
w 


= 


w 

£ 


2E. = 


-0.000174 


N* ' 

V 


= 


• > 


II 


-0.000176 


Z* ' 

q 


= 


-M* ' 

w 


= 0. 


000174 


M* ' 

q 




-Z* ' 
w 


(¥1 


) 

= -0.000091 


N- ' 
r 


= 


Y* ' 

V 


i¥f 


= -0.000093 


K*’ 

V 


= 


Y* *2 

V 


c£. _ 


-0.000020 
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(2.32) 

(2.33) 

(2.34) 

(2.35) 

(2.36) 

(2.37) 

(2.38) 



I 



y.» = N*' = 0.000176 

r V 



(2.39) 



Contributions due to the deflections of the fins result in the 
following coef f icient Sj for which it suffices to calculate those for one 
fin only. 

' ■ 



a = = 3.394 

C, = 3.95 

La 

= 46.8 ft^ 



YI = 0.005706 

0 



Z 



6 




-0.007656 



2 

where = 3.6 and A^ = 68 . 9 f t . 

Y' |Z" 1 

K' = = 0.000269 

0 Z 



(2.40) 



(2.41) 



where Z' ' is the distance for the moveable portion, -8.46 ft. 
cp 



M' = = -004019 

0 36 



(2.42) 
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i 



where = 28.8 ft^. 






„2 ^ Z 



N. = -0.001680 

0 



(2. A3) 



When the fins are deflected, drag forces are induced. The drag 
force of a body is given by 

1 o .2 

X = — Cj^pAu 



From Comstock [8] , 



C 



D 



DO 




irae 



where 



a = AR^ = 1.182 

e = 0.9 (Oswald efficiency factor) 



Since 



C 



L 




1.656 



then 



"D 



= CjjQ + 0.81A66' 
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Substituting this into the X equation, 



X = 



+ 0-8146 



^)pA^ 



and, taking the second derivative with respect to 6 one gets 



= 0.8146PA^U^ 



Next , non dimens ionalizing. 





-0.004053 



(2.44) 



2.3.3 Sail Coefficients 

Figure 5 shows the dimensions of sail used. Using the method 

described in Reference 8, to evaluate C., for the sail one calculates 

La 

the following. 



A = 0 radian for sail 



AR = a 
s 



( 2 ) ( 10 ) 

16 



= 1.25 



A = 160 ft 



27T 



La 






= 2.42 
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SCALE: 1/4 in. = 1 ft 




h< 8 ft H 

58 ft VW-*^ F.P. 




SURFACE AREA = 



27T Vl/2(a^+W* h 
397.4 ft2 



Figure 5. Sail geometry. 
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Y' 

V 



-0.011934 



(2.45) 



The center of pressure is assumed to be located at the 1/4 chord length 
from the forward edge. 



X = 28 ft 

cp 



-Y 'x 

N' = - = 0.001856 (2.46) 

V 5, 



The contribution to roll is 



K’ 

V 



Y’z 
V cp 

I 



where z is the distance from the centerline of the vehicle to the 
cp 

vertical center of pressure, and z = z’ = -15 ft for this submersible. 

cp cp 

Then, 



K’ = 0.000155 (2.47) 

V 

Y’ = N’ = 0.001856 (2.48) 

r V 

Y ’x 

N’ = = 0.000289 (2.49) 

r £ 

K' = Y' (-^) = -0.000083 (2.50) 

p V \ I / 
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The acceleration derivatives, according to Abkowitz can be de- 
rived analogously as for the stern control surfaces. In this case, the 
span will be twice the height of the sail and the resulting derivative 
will be divided in half and then nondimensionalized . In the vertical 
plane one gets 



- 



y. » 

V 



2 2 
ITpS C 



0.54 

1 + - + - 
c s > 






Y* ' 

V 



= -0.000886 



where 



s = 20 ft 

c = 16 ft 



Hence 



-Y* 'x 

N-' = = 0.000138 

V I 



K- ' = 

V 



Y* 'z 



= 0.000074 



(¥f 



N‘ ' = Y 

r V 



= - 0.000021 



Y-' = N*’ = -0.000138 

r V 






K- ' = Y 

P 



= -0.000006 



(2.51) 



(2.52) 

(2.53) 

(2.54) 

(2.55) 

(2.56) 
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2.3.4 Sailplane Coefficients 



Figure 6 shows the geometry of the sailplane. The center of 
pressure is assumed to be 1/2 mean chord length forward of the trailing 
edge. For this vehicle. 



X = 24 ft 

cp 



The following definitions are used in calculating A and C 



La 



b = tip of fin to hull centerline 

= 11 ft 

d = maximum diameter of hull 

= 20 ft 

a = aspect ratio = 2.67 for a single plane 

A = 31.5 ft^ 

^La = 3-59 

-C A 

Z' = — ^ = -0.003490 

-Z’x 

M' = — = 0.000465 

w a 

There is no contribution to roll in the horizontal plane. 



(2.57) 



(2.58) 
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SCALE: 1/4 in. = 1 ft 




K2ft-H 




PROJECTED AREA OF 
SINGLE FIN = 31.5 ft2 



Figure 6. Sailplane geometry. 
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0.000465 



(2.59) 



Z' 

q 



w cp 

I 



M' 

q 



Z'x 

q cp 

I 



0.000097 



(2.60) 



The contributions due to the deflection of the fins are calcu- 
lated in a similar fashion as before. 



7 ' 

6b 



Ct a 

La 



= -0.005429 



where z'' for the sailplane is 16 ft. 
cp 



M' = 

6 



Z 1 1 X I 
6 cp 



= 0.000724 



(2.61) 



(2.62) 



As with the stern control surfaces, when the fins are deflected, 
drag forces are induced. 

a = AR = = 1.2727 

msp - 

e = 0.9 (Oswald efficiency factor) 



Since 



C 



L 




3.206 
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then 




^DO 




TTae 



= CjjQ + 0.450162 



Substituting this into the expression for the drag force of a 
body X = *2 Cj^pAu gives 



X = 



+ 0.45016^^pA___u^ 



msp 



Taking the second derivative with respect to 6, 



X.. = 0,4501pA u 

66 msp 



Nondimen sionali zing. 



-X 



''66 



66 



p£' 



■ 2 \ 

N 



= -0.001750 



(2.63) 



The acceleration derivatives for the sailplanes are calculated as 
follows, with c = 5.25 ft and s = 14 ft. In the horizontal plane, 



1 - 



0.54 



Z* ' 
w 



ttp(14)^(5.5)^ 

4V(14)^ + (5.5)^ 



1 + 



14 5.25' 



5.25 14 
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0.000092 



(2.64) 



Z- ’ = 

w 



M* ’ 
w 


_ ^ 

1 


= 0.000012 


( 2 . 65 ) 


Z* ’ 

q 


= M- ’ 
w 


= 0.000012 


( 2 . 66 ) 


M- ’ 

q 


Z- ’ 
w 


= - 0.000002 


( 2 . 67 ) 



The hydrodynamic data to be utilized by the SUBRUN computer program, 
that will generate the linear dynamics, requires that the coefficients 
be input in a dimensionalized form rather than a nondimensionalized 
form. Appendix D contains a listing of dimensionalized hydrodynamic 
data and dimensionalizing factors for this purpose. 

2.4 Computer Program Description 

The computer program, SUBMODEL, used for this thesis was devel- 
oped at Draper Laboratory [9]. The program has the capability to inte- 
grate the nonlinear equations of motion of a submarine; search for a 
local equilibrium point in the nonlinear equations of motion; calculate 
the linearized dynamics of the linearized equations about a nominal 
point; integrate the linearized equations of motion; simulate the com- 
pensated vehicle using either the nonlinear or linear equations. 

As previously stated, the nonlinear equations take the form of 
Equation (2.1) 



Ex = f_(x, u) 
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where 



X = 10 X 1 state vector with states u, v, w, p, q, r, phi, 

theta, psi and z, the last four describing the vehicle 
attitude with respect to an inertial reference frame. 

u = 4 X 1 control vector with DB, Dl, D2, D3 as members 

^ = 10 X 1 vector that is a nonlinear function of the states 
and the controls 

^ = 10 X 10 matrix 

The nonlinear equations were integrated with the RPS propulsion 



model, as described in Appendix E, using as input a velocity in feet/ 
second (corresponding to 21 knots) with initial values of zero for all 
states, with the exception of Dl deflected to 2 degrees, D2 deflected to 
0.5 degree and D3 deflected to -0.4744 degrees. The assumption is made 
in the propulsion data that wake and thrust deduction factors vary 
slightly and are therefore considered constant. The nominal point is 
then determined by integrating the nonlinear equations with the eta 
propulsion model using as inputs for the states, the final value of the 
states from the RPS integration. 



Once the nominal point has been determined, the linearized dynam- 



ics can be calculated by linearizing about that nominal point (see Ap- 
pendix E for a description of Nominal Point). Recall that the linear- 
ized equations are of the form 



EAx = A’Ax + B’Au 



( 2 . 68 ) 



or 



equivalently, the conventional representation form 



X 




(2.69) 
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On can easily see that A in state space form equates to E and ^ 
equates to ^ . It is E ^A’ and ^ that define the linear system 

dynamics, on which the control design is based. 

To validate the resulting linear model, one perturbs the states 
of the nominal point and compares the resulting responses of the per- 
turbed linear and nonlinear models. Perturbation from the nominal point 
on the order of 10% is adequate for checking validity of the linear 
model behavior in comparison to the nonlinear model. Figure 7 shows the 
graphical responses of the two models for comparison. 

Once the controller has been designed^ the effectiveness of its 
operation can be tested by using it on the previously derived linear and 
nonlinear models and analyzing their responses to step and ramp inputs. 

A description of the controller development is given in Chapter 3 with 
the results of the subsequent testing in Chapter 4. 
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Figure 7. Linear and nonlinear submarine model responses. 
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Figure 7. (Cont.) 
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Figure 7. (Cont.) 



44 



NON-LINEAR RESPONSE (CONTD) 




TIME 




TIME 




time 




TIME 




Figure 7. 
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CHAPTER 3 



COMPENSATOR DESIGN 



3.1 Overview of the LQG/LTR Procedure [10] 

In this chapter, a controller will be designed using the LQG/LTR 
design procedure. The singular value loop shaping approach will be used 
to mold the singular values of the system open-loop transfer function to 
meet the specifications of performance and robustness to plant uncer- 
tainty and modeling errors. In order to meet the performance specifi- 
cations of small steady state errorsj one desires high DC gains. To 
minimize the effect of wave encounter and modeling errors at high fre- 
quencies, we desire crossover at about 0.1 radians per second with a 
large roll-off after crossover. This will ensure attenuation of high 
frequency modeling errors and sensor noise, as well as wave effects 
which typically occur in the 0.2 to 2 radian range. 

The LQG procedure is a mathematical method of designing a robust 
multivariable controller for a state space system. The designer, how- 
ever, has the flexibility to pick the parameters of the design motivated 
from a frequency domain setting. An iterative design process is then 
carried out, during which the designer is able to match desired fre- 
quency domain characteristics, captured by the singular values of the 
multivariable system, which are equivalent to the SISO Bode plots^ by 
using time-domain mathematics as derived for the solution of the origi- 
nal optimal control problem. In this way, a systematic procedure can be 
carried out, with the mathematical solution kept simple and in closed 
form for even a complex MIMO design. 
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The design is based on the feedback loop as represented in Fig- 
ure 8. The plant ^(s) has an input _u, and output a disturbance d as 

seen at the output, and measurement noise vector n. The compensator is 
designated as K(s) . 



d 




Figure 8. Generalized feedback loop. 

The plant G is the physical system to be controlled and can 
therefore behave in a complex manner, as implied by the term nonlinear, 
infinite dimensional and time-varying which are used for its modeling. 

To simplify the situation, we desire G to be linear, finite and time- 
invariant (as results from the linearization procedure described in the 
preceding chapter), with a transfer function ^(s). 

The design objective will be to provide good command-following, 
good disturbance rejection, and small responses to sensor noise, all 
subject to the constraints imposed by modeling errors. Hence we need 
large loop gains in the frequency ranges where commands and disturbances 
are large and small loop gains in the frequency ranges where modeling 
errors are large, as well as a well-behaved crossover. These frequency 
domain requirements can be captured as shown in Figure 9, where the 
maximum and minimum singular values of the loop transfer matrix envelope 
the MIMO system performance, and can be directly interpreted as SISO 
Bode plots. For a good design, it becomes necessary to use techniques 
which directly take plant descriptions and performance specifications 
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Figure 9. Frequency domain requirements. 



into consideration and systematically produce a controller which satis- 
fies the singular value requirements of Figure 9. The multivariable 
design technique, LQG with Loop Transfer Recovery, is able to perform 
such a direct singular value loop-shaping. 

As previously stated, the LQG/LTR methodology was developed to 
perform control system design on the basis of time-domain optimization. 
The method uses the following state space description of the plant , G: 

X = Ax + Bu + 5 (3.1) 



Z “ ^ H 



(3.2) 



where x is an n-dimensional state vector, u and y_ are m-dimensional 
control inputs and outputs, and and t\_ are white noise processes. The 
A, and matrices satisfy 



G(s) = C^(s)B 



(3.3) 
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where ^(s) = (S^ - A) . The feedback design task is to find a control 
that satisfies the specifications. For this purpose, the math- 
ematics developed for the solution of the following optimization problem 
are employed^ although the design is geared towards loop— shaping and 
clearly cannot claim optimality as will be seen in the sequel. 

T T T 

Min J = J [x (t)^(t) + u (t)^(t)] dt (3.4) 

0 

where the R matrix is the control weighting matrix and Q is the state 
weighting matrix. The solution of this optimization problem is given by 

u(t) = -Gx(t) (3.5) 

where ^ is a full-state linear quadratic regulator (LQR) gain as defined 
in the Control Algebraic Riccati equation (CARE). However, in the LQG/ 
LTR procedure, the solution is carried out in two stages. First, the 
filter gain matrix, H, is found by solving the Filter Algebraic Riccati 
equation (FARE) . 

The objective in this first step is to shape the singular values 
of the Kalman filter (filter gain) matrix according to the specifica- 
tions, as suggested in Figure 9, and then recovering the same loop 
transfer matrix shape for the singular values of the overall system loop 
transfer matrix ^(s), which includes the plant, any augmentation, and 
the controller gain matrix G, which is obtained as a second stage solu- 
tion, by solving the CARE, for a given filter design. This compensator 
structure is referred to as a Model Based Compensator (MBC), because it 
makes explicit use of the plant linear model (i.e. A, £ matrices). 

At the end of the second stage, the overall system matrix is obtained as 

T(s) = Gp(s)K(s) (3.6) 
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where ^(s) — — A + BG + HC ) H, Once the shape meets the design 

specifications, we are further guaranteed that the poles of [A - HC] and 
[A - BG ] will lie in the left half plane, i.e., no instabilities will be 
introduced into the nominal system, provided [A, B] are stabilizable and 
[A, detectable. 

3.2 Modal Analysis 

Since this thesis emphasizes primarily the influence of the MDS 
configuration, the decision was made not to address the effect of the 
sailplanes. Additionally, the states psi and z are of no direct inter- 
est because of their non-effect on the other (body) states. The origi- 
nal plant A matrix contained zeroes in the columns designated for psi 
and z. By ignoring psi and z, this reduces the plant matrix to an 
(8 X 8) matrix, the control matrix to an (8x3), and the output matrix 
to a (3 X 8) . It then becomes necessary to determine which states to 
use for the purposes of control. 

The MDS arrangement has the unique potential to control roll, (}) , 
(and pitch, 0) via D2 and D3. Finally, for the third variable, turning 
velocity, r, (yaw rate), could be controlled with Dl. 

At the outset, the corresponding poles and zeroes of the plant 
were determined. The variable combination (|), 0, and r provided no non- 
minimum phase zeroes and no unstable poles. Table 1 lists both poles 
and zeroes for the open-loop plant. 

A modal analysis was then conducted to ascertain the relative 
controllability and observability of the system. To examine the con- 
trollability of a system, one is interested in the ability of (a) parti- 
cular input (s), in this case, Dl, D2 and/or D3 , to control the (entire) 
state of a system. The observability of a system tests the possibility 
of observing the complete state of the system through a particular out- 
put, i.e., r, (p or 0. Using the augmented system A, ^ and matrices 
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Table 1. 



Augmented Plant Poles 

-2.029353579E-02 

-4.377357711E-02 

-1.803910515E-01 

-1.843174783E-01 

-6.227811352E-01 

-6.583723063E-01 

-5.085034079E-01 ± 5.224757797E-01 
0.000 (Multiple of 3) 

-l.OOOE-01 (Multiple of 3) 

Augmented Plant Zeroes 

-4.421289192E-02 
-1. 98463495 lE-01 
-2.583928560E-01 

(augmentation to be described in Section 3.3), the corresponding eigen- 
values and eigenvectors were calculated and a modal analysis performed 
in order to study the controllability and observability properties of 
the system. The eigenvalues give the response time of a system's modes 
while the eigenvectors represent the relative contributions of each 
state to a particular mode, represented by the corresponding (normal- 
ized) eigenvector. For example, the plant can be described in state 
space from by. 



X = Ax + Bu (3.7) 

^ ^ (3.8) 



51 



( 

I 



I 



! 



If we let X Tz , for any nonsingular we can define a new system. 
However, one must remember that £ does not describe the physical state 
space. ^ can be chosen to be the matrix of eigenvectors. Using this 
change of variable transformation, one gets that x = Tz. So if one 
makes this change of variable substitution into the state space equa- 
tions (3.7) and (3.8), one is left with 

X = Tz = ATz + Bu (3.9) 

Z = (3.10) 

To complete the transformation, T ^ is multiplied through (3.9) and 
(3. 10) thus yielding 



£ = T ^ATz + T ^Bu (3.11) 

The net effect of the above procedure is the diagonalization of the 
state equation where the elements on the diagonal correspond to the 
eigenvalues of the system, if they are distinct. Rewriting the modal- 
domain state and output equations for MIMO system, one gets 

_z = Az + B \ (3.12) 

Z = c\ (3.13) 

* _X * 

where ^ = T B, _C = CT, and ^ is a diagonal matrix of eigenvalues. 

For a system to be controllable by an input _u or observable by output 

:k "k 

the column matrix B or the row matrix ^ must not contain zeroes in a 
SISO setting. In a MIMO setting, if zeroes appear, then there must be 

at least one nonzero entry in the B matrix and one nonzero entry in the 



52 



( 



* 

matrix for a particular state to be controllable and observable by at 
least one corresponding input or output. 

Another method could be used to determine the degree of control- 
lability and observability for a MIMO system and that is the rank meth- 
od. If one forms the matrices 

[B, AB, A^B, a""^B] 

and 





and their rank is equal to the order of the system, then the system is 
both controllable and observable. If the rank is less than the order of 
the system, then that difference is the number of uncontrollable and 
unobservable states [11]. Appendix F contains the eigenvalues (exclud- 
ing augmentation) and the corresponding eigenvector matrix, TEMP. The 
matrix is the diagonalized modal domain matrix with the eigenvalues 

•k k 

located along the main diagonal. The new B and ^ matrices are also 
included and one can observe that ^ zeroes exist in either matrix. 

Thus one can conclude that control inputs Dl, D2 , and D3 can each con- 
trol all the states of the plant while the entire state vector can be 
observed by any of the r, (p, and 0 outputs. 
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3.3 



System Augmentation 

Before setting off to design a compensator for the control sys- 
tem, one must examine the behavior of the plant singular values. Fig- 
ure 11 shows the open loop singular values with the system augmented 
with an integrator for each input. In order to further improve the 
model, one must take into account the actuator dynamics. Actuators are 
normally modeled as second order systems, however, when one takes into 
account the high damping in a real actuator, one can then approximate 
the dynamics of the actuator as a first order system. Pictorially this 
situation looks like Figure 10. 

The state-space equations for this open loop representation are: 



• 

X 


= Ax + B u 


(3.14) 


—p 


— p-p — p-p 




y 


= C X 


(3.15) 




-p-p 




X 


= Ax + B u 


(3.16) 


—a 


—a— a — ar-a 






= C X 


(3.17) 


—a— a 




X^ 


= A^x,. + Bt-Ut 


(3.18) 


-I 


—I— I —I— I 






II 


(3.19) 
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Combining the above and putting into matrix form, one gets 



X 




B C 
A 

— a 
0 










0 




0 ] 



X 

—a 






The modeling of each actuator to a first order approximation will 
be of the form 



G 

a 



1 

TS + 1 




X 

a 

u 



( 3 - 20 ) 



Therefore 



Sx 

a 




u 

T 



( 3 . 21 ) 



In state space form, this will be represented as 



X 



a 




( 3 . 22 ) 
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For a system with three control surfaces 



in matrix form, we get 



■_ 

T 


0 


0 




T 


0 


0 


0 


T 


0 


B 

—a 


0 


T 


0 


( 

O 


0 


_ ^ 
T_ 




0 


0 


T . 






1 0 0 
0 1 0 
0 0 1 



The addition of three integrators is needed to provide high 
enough gains at low frequencies for good command following and disturb- 
ance rejection. 



G 



I 





u 



( 3 . 23 ) 



Thus for each control variable x^s = u or x^ = u. Then for a system 
with three control inputs, the associated matrices will look like 





"0 


0 


o " 




"1 


0 


0 " 




0 


0 


0 


ltd 

M 

II 


0 


1 


0 




.0 


0 


0 _ 




_0 


0 


1 . 







0 

1 

0 



0 

0 

1 
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The. augmented system in state-space form is consequently given by 





















— 1 


00 

X 

00 

<\ 


— I>-a 




(8 


0 


3) 






0 






(8 X 3) (3 X 


3) 


X 




(8 


X 


3) 




0 


0 


J. 


0 


0 












T 




X 










0 


0 -i 


0 


0 




0 


X + 




0 




(3 X 8) 


X 






X 






(3 


X 


3) 




0 0 


_ 


0 


0 
















T 






X 
















0 


0 


0 




1 


0 


0 


0 


0 




0 


0 


0 




0 


1 


0 


(3 X 3) 


(3 X 3) 




0 


0 


0 




0 


0 


1 



(3.24) 



Z 



"0 0 
0 0 
0 0 



0 0 
0 0 
0 0 



0 1 
0 0 
0 0 



0 0 
1 0 
0 1 



0 

(3 X 3) 



0 

(3 X 3) 



Note that the matrix was not included in the output matrix because 
the desired outputs are r, c(), and 0. 



The time constant, t, is defined as — ~r [12]. 



Recall that the 



projected crossover frequency is about 0.1 radian per second so as to 
avoid the natural frequency of ocean waves of 0.14 radians per second. 

A typical value for the natural frequency and damping ratio for an actu- 
ator is 3.14 radians per second and 5 = 0.9. The resultant time con- 
stant for the actuator then becomes 0.354 sec. This is clearly above 
the desired bandwidth. It was consequently decided to ignore the actu- 
ators and instead use lag compensation at 0.1 radians per second to 
improve roll-off near crossover. This will have the effect of making 
the vehicle less susceptible to modeling errors and high frequency 
noise. Therefore, the system is augmented with three lag compensators 
whose dynamics are as shown in Figure 12. 
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Figure 10. Augmentation representation. 




Figure 11. Singular values of open loop plant. 





0.1 






s + 0.1 





Figure 12. First order lag compensator. 
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3.4 



Controller Design 

The method employed to design the controller uses the Model-Based 
Compensator approach [13]. Two very important points must be kept in 
mind. First, the controller must be able to maintain the stability of 
the closed loop system and, secondly, be able to achieve desirable sin- 
gular value loop shapes, in the frequency domain which translate to 
performance specifications. This can be achieved by appropriately shap- 
ing the system singular values in a manner analogous to the single in- 
put-single output Bode plot shaping. Pictorially^ the controller and 
system look like Figure 13. Figure 14 illustrates the state-space de- 
scription of the Model-Based Compensator, K(s) , and the open-loop plant, 
G (s). It should be noted that the A, B, and C matrices that appear in 

-p 

the plant also appear in K(s) in a similar fashion and thus the term 
"model-based . " 



d(s) 




Figure 13. Feedback structure of a MIMO control system. 
The dynamics of the open loop plant can be written as 



x(t) 


= Ax(t) + B_u(t) + Ld(t) 


(3.25) 


z(t) 


= Cx(t) 


(3.26) 



The open- loop transfer matrix is then 

= C(S^ - A)“^B 
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Gp(s) 



(3.27) 



so that 



X(s) = Gp(s)u(s) (3.28) 

The Model Based Compensator dynamics utilizes the state vector z(t) 
which has the same dimension as the plant state x(t) . The dynamics of 
the MBC can be written as 



_z(t) 


= ^(t) + ^(t) + Hv (t) 


(3.29) 


v(t) 


= -e(t) - Cz(t) = x(t) - C£(t) - £(t) 


(3.30) 


u(t) 


= -Gz(t) 


(3.31) 



In the above, the augmented dynamics of the plant as described by 
Eq. (3.24) are used. Combining the above three expressions we get 

£(t) = [A - ^ - HC]^(t) - He(t) (3.32) 

Thus the input to the controller is the error, ^(t) and the output is 
^(t) . The input-to-output relation can be written as 

u(s) = K(s)e(s) (3.33) 

where K(s) = G(s^ ~ A uses an appropriate change 

of variable, one can determine the stabilizing nature of the gain ma- 
trices G^ and H. The appropriate variable change is 

w(t) = 2 l(t) - £(t) (3.34) 

It follows that w(t) = x.(t) - _z(t) . After making the appropriate sub- 

• • 

stitution for x(t) and z(t), one is left with 
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w(t) 



(3.35) 



= [A - _^]w(t) + ^(t) - l^(t) 

This resulted in the decoupling of the dynamics of w(t) from that of 
x(t) . Then x(t) can be written as 

x(t) = [A - ^]x(t) - ^(t) + ^(t) (3.36) 

The dynamics of the closed-loop system in matrix form is then given by 




The eigenvalues of the closed loop matrices [A - and [A - K] con- 
stitute the system poles. Therefore, the det (X^ - A + BG ) and 
det (X^ - AH- HC ) provide the needed insight into the stability of the 
closed-loop system. In order for the closed— loop system to be stable 

ReX^[A - ^] <0; i = 1, 2, ..., n (3.37) 



and 



ReX^[A - HC ] < 0; i = n+l,nH-2, ...,2n 

If we are given matrices A, ^ and ^ with [A, B] stabilizable and [A, C] 
detectable, then matrices ^ and H exist that stabilize the open loop 
system, i.e. [A - HC ] and [A - are stable. The methodology used to 
calculate ^ and H is outlined in sections 3.5 and 3.6 [13, 14]. This 
provides in fact the underlying mathematical basis of the LQG/LTR ob- 
tained solutions. 
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3.5 



The Filter Gain and Control Gain Matrices 



It can be seen from Figure 14 that the filter and control gain 
matrices are the heart of the MBC. In essence the Kalman filter gain 
matrix, H, is the input matrix into the MBC taking the error input 
through the MBC dynamics, which in turn passes it through the integrator 
having first multiplied the outputs by the appropriate gain matrix G. 
This then generates the appropriate control input for the plant. 

The G and II matrices serve a very useful purpose. Recall from 
Section 3.1 that the desired plant output is modeled via singular val- 
ues. Setting the shape of the desired plant open loop singular values 
is important for good command following, disturbance and noise rejec- 
tion. Once the designer has achieved a shape that meets the open loop 
specifications, it is important to realize them next in a closed loop 
configuration. The feedback configuration has the advantages of reduced 
sensitivity to external disturbances, as a result of a course change — 
for example, depth change, as well as meeting the specifications of good 
command following and other noise and modeling uncertainty reduction. 
Through the clever iterative selection of the H and ^ matrices via the 
LQG/LTR procedure this can be attained. Section 3.6 will detail how 
this selection is made. 

3.6 Loop Shaping 

This section will outline the procedure that was followed in 
designing the LQG-based controller with Loop Transfer Recovery [14]. A 
plot of the open loop singular values of - A) ^ is given in Fig- 

ure 11. Since there are three controls and three outputs, there are 
three lines of singular values. As one can see, there is significant 
separation between the maximum and minimum singular values. It is de- 
sirable that all three behave as close as possible in a tight controller 
design. 
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The first step is to shape the open loop plant singular values so 
they meet the desired specifications of good command followings stabil- 
ity and performance, as well as noise and disturbance rejection. This 
is done by tying the maximum, middle, and minimum singular values to- 
gether at high and low frequencies. To achieve this, one must form the 
matrix 



G^olCs) = C(SI - A)“^L (3.38) 

where L is an arbitrary matrix that can be manipulated by the designer 
to ultimately give the desired loop shaping. For system at hand, we 
form. 



(SI - A) 



(SI - A*) 

0 SI 



where 



A 

— X 



A 

(8 X 8) 
0 

(3 X 3) 



B C 
— p — a 

(8 X 3) 




0 



0 




0 0 



0 

0 

I 

T 
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.1 



and 



1* 



0 

(8 X 3) 



Note, the A matrix of (S^ - A) came from Section 3.3. The inverse of 
(SI - A) is 



(SI - A) 



-1 



(SI - A^) 



-1 

1" (SI - A^) 
1 
s’ 



-1 



The L matrix can be broken up into two components and where 

For low frequency shaping, SI is 






-(8 X 3 ) ^2 



-(3 X 3 )' 



approximately zero. With that approximation and forming G (s) , one 

rUL 

gets 





-1 1 






-1 -A* 




C(SI - A)’^L = ^ 0] 


-A. 




I 

0 z 


—1 




L ~ s J 





(3.39) 



where ^ is a (3 x H) matrix with identity in positions (1,6), (2,7), 
and (3, 8) corresponding to the states to be controlled r, 4), and 0. 
Expanding (3.39) yields 
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C(SI - A) = 



„ .-U - 11 ^ -’"^2 

-1 ■ S 



(3.40) 



But for small s, low frequency, the second term in (3.40) dominates the 

- 1 1 lA* — 

first term. So ^(SI_ - A) L can be approximated by r . Then 

to tie the singular values together at low frequencies we also have 
1 

Gfol(s) = -g- Therefore 



k2 ' 



( 3 . 41 ) 



For high frequencies, one can approximate (SI_ - A) by S_I. Using 
the same procedure as above. 



C(SI - A) 






0 ] 



(SI) 



-1 



- (SI) 



-1 



^2 



(3.42) 



Expansion of (3.42) yields 



C(SI - A) 



- 11-1 - ll -*^2 



(3.43) 



As S gets large, high frequencies, the second term in (3.43) will be 
dominated by the first term. So will be equal to 



h - (£u£ii)''£u - £L(^ii£u)'‘ 



(3.44) 
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The combination of and constitutes the matrix of the 

overall system. The matrix provides loop shaping and is totally under 
the designer’s control as a free design parameter. Figure 15 shows the 
result of the initial open loop shaping. This, in turn, provides the 
desirable system profile which the closed loop compensated system should 
match. The resulting L matrix can be found in Appendix G. In order for 
the crossover frequency to be less than the critical crossover frequency 
of 0.1 radians per second, (natural sea spectrum is at 0.14 radians per 
second), it was necessary to multiply the L matrix by a scalar of 0.05. 







Figure 15. Singular values of 

The next step in the controller design is to determine the filter 
gain matrix, H. To do this one solves the Filter Algebraic Riccati 
Equation (FARE). 
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0 



1 T -1 
-J- ZC £ CZ 



(3.45) 



= ZA + A^Z + 5 - 



where ^ is an arbitrary (m x m) symmetric and positive definite matrix 

(equal to the identity here) and normally the process noise covari- 

T 

ance matrix is now replaced by . The solution to the FARE is the 
(n X n)^ matrix. The filter gain matrix is given by, 

1 T -1 

H = -^5 (3.46) 

If one plots the singular values of the Kalman filter, 

= C(SI - A)"^H (3.47) 

then they should pattern the profile of the singular values of G (jco) . 

— r UL 

Figure 16 shows the similar profile. The DC gains, crossover frequency, 
and band crossover bandwidth are approximately the same. If it becomes 
necessary to adjust the DC gain and crossover, this can be accomplished 
by selecting another free parameter y , in Eqs. (3.45) and (3.46). 

The final step in the controller design is to solve the Control 
Algebraic Riccati Equation (CARE) to obtain the control gain matrix, G. 
The CARE looks like 



0 = - a\ - £Q + KBR~^b\ (3.48) 

where K is an (n x n) solution. The R matrix is an (m x m) positive 
definite matrix (identity for our purposes) , and the Q matrix is given 
by 



Q 



T 

= C C 



(3.49) 
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SINGULAR VALUES (LOG MAGNITUDE) 

- 3.75 - 2.50 - 1.25 0.00 1.25 2.50 3.75 




Figure 16. Singular values of G^(jo)). 

such that [A, is a detectable pair. The control gain matrix G is 
given by 



-1 T 

G = R B K 



(3.50) 



The singular values of the resultant loop transfer matrix 
(o^[T^(jai) ] ) are now calculated where 



T(s) = G(s)Kj^q^Cs) (3.51a) 

or 

T(s) = Gp(s)K(s) (3.51b) 
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£(s) is the plant transfer matrix defined by 



^(s) = ^p(s)Ga(s) 



(3.52) 



as illustrated in Figure' 7. 
by: 



In the frequency domain ^qq(s) is given 



^Qg(s) = G(S^ - A + ^ + HC) (3.53) 

Figure 17 is a plot of the singular values of the compensated system 
T^(s) . A value of 1000 was used for q so that the three individual con- 
trol singular values coincide. It is apparent by comparing a^(^(jo))) 
with cr.(^ (j(jj)) that the DC gains are equal and crossover occurs at the 
same frequency. The roll-off past crossover is desired as illustrated 
in Figure 9. 

A check was made to determine if the controller, K(s) introduced 
any instabilities into the system. Table 2 is a listing of poles and 
zeroes which indeed corroborates the fact that the design has resulted 
in a stable system. Appendix H contains the A, ]B, H and 

(A “ ” HC ) design matrices that will also be used in the simulation 

as described in Section 4.1. 
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SINGULAR VALUES (LOG MAGNITUDE) 

-12.0 -B.o -4.0 0.0 4.0 8.0 12.0 




Figure 17. Singular values of T(ju). 
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Table 2. Poles and zeroes of K(s) . 



Unique Poles 
Number Real 

1 -4.421317785E-02 

2 -1. 984642 130E-01 

3 -2.515003179E-01 

5 -2.283572447E-01 2 

7 -3. 121226887E-01 7 

8 -8.357975180E-01 

10 -8.470467465E-01 3 

12 -4.856380807E-01 1 

14 -1.2022 1261 9E+00 5 

Unique Zeroes 

Number 

1 -2.418544925E-02 

2 -4.342940086E-02 

3 -9. 122782906E-02 

4 -1.062659914E-01 

6 -1.471524383E-01 -7 

7 -1.852788869E-01 

8 -5.245929103E-01 

9 -6.230977444E-01 

11 -5. 179401036E-01 4 



Imag 



689624484E-01 

200466481E-01 

140172344E-01 
221 1596 15E+00 
654424752E-01 



231643907E-02 



.323091793E-01 
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CHAPTER 4 



COMPENSATOR PERFORMANCE 



4. 1 Model Simulation 

The model simulation used to test the controller involves the 
integration of 



£(t) = ( A - ^ - ^)^(t) - He(t) 

where ^(t) is the error vector defined by 

e(t) = _r(t) - x(t) 

To close the loop, one can either select the linear model or nonlinear 
model as developed in Chapter 2. It is expected that the linear model 
will give good results with the nonlinear model giving acceptable re- 
sults because the controller was designed from the linear model origin- 
ally. 



4 . 2 Simulation Results 

Now that the controller has been designed, it is necessary to 
test its performance. Recall that the motivation behind this thesis was 
twofold: (i) design a controller from a linear model and (ii) test it 

on a nonlinear model, controlling r, (j), and 0 in a turn. 
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The first test simulation run purpose was to check if the con- 
troller would be able to achieve zero steady state error with a command 
reference input of -0. 368482E-02 rad/s, the nominal r value, and zero 

for both (j) and 9, keeping in mind that the open loop roll for the same 

-2 

condition was 1.5 degrees and 9 4.0 x 10 degrees. With the closed 
loop simulation, both (j> and 9 indeed went to zero. 

A second simulation run was made to command a condition that was 
other than the nominal r. A turning velocity equivalent to a 5 degree 
nominal rudder deflection was commanded with zero again for (() and 9. 

The open loop (}) and 9 for this condition were 5.5 degrees and 2.5 de- 
grees respectively with the controller, r went to the commanded value of 
-0.1160E-02 rad/s^ and (}) and 9 went to zero. With both simulation runs, 
the controller achieved zero steady state error for the three command 
inputs. Figure 18 contains the open loop r, (J) and 9 response runs and 
the closed loop simulation response with command input of -0 . 368482E-02 
rad/s, 0®, and 0® for r, (}) and 9 respectively. Figure 19 likewise shows 
the open and closed loop responses for r = -0.1160E-01 rad/s and zero 
for (j) and 9 . 

One additional benefit of the controller can be seen with depth, 
z. Without the controller, both open loop simulations exhibited contin- 
uous depth change throughout the run. With the controller^ the vehicle 
sank a small amount but steadied as the run continued. Therefore the 
controller was able to stabilize depth loss even though depth was not a 
control variable. 
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Figure 18. Open and closed loop simulation responses with 
command input . 

r = -0.368482E-02 rad/s 

<f) = 0® 

0 = 0 ° 
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Figure 18. (Cont.) 
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CLOSED LOOP SIMULATION RESPONSE (CONTD) 





TIME 




TIME 

Figure 18. (Cent.) 
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OPEN LOOP RESPONSE 




TIME 






Figure 19. Open and closed loop simulation responses with 
command input . 

r = -0.1160E-01 rad/s 

(j) = O*’ 

0 = 0 ° 
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Figure 19. (Cont.) 
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Figure 19. (Cent.) 
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CHAPTER 5 



CONCLUSIONS AND RECOMMENDATIONS 



The LQG/LTR method can be successfully used to design a control- 
ler from a linear model. Additionally the resulting controller is able 
to perform acceptably with zero steady state error on a nonlinear model, 
which more closely approximates the real world system. 

The MDS configuration has demonstrated improved ability over the 
conventional cruciform stem in turning situations. With the ability of 
the vehicle to align its local reference axis to that of the global 
system, one is able to obtain true motions of surge, sway^ heave, roll, 
pitch, and yaw. The advantage being improved crew comfort and quicker 
computer solutions to fire control problems. 

It is recommended that further thesis work on this model to in- 
clude depth control with the additional control surface of sailplanes. 

A potpourri of operating points spanning a speed range from 5 to 21 
knots with rudder deflections of 5 to 30 degrees will give on a feeling 
for the capabilities of the MDS configured submersible. The addition of 
vortex shedding to the model will also improve the characteristics of 
the model to real world submersibles . Finally, it would be desirable to 
design a controller that could negotiate ordered changes in depth. With 
gain scheduling, the model could clearly be used as a training aid in a 
simulation laboratory . 
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APPENDIX A 



MASS PROPERTIES 



Mass 

75,483 slugs 
Ixx 

3,019,300 slugs-ft^ 
Ixy 

0 . 

XG 

0 . 

XB 

0 . 

Xstern 
-89.75 ft 



Weight 

2,428,300 Ib^ 

lyy 

123,790,000 slugs-ft^ 
Ixz 

0 . 

YG 

0 . 

YB 

0 . 

S length 
180 ft 



Buoyancy 

2,428,300 Ib^ 

Izz 

123,790,000 slugs-ft^ 
lyz 

0 . 

ZG 

1,0 ft 
ZB 

0 . 
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APPENDIX B 



NONLINEAR EQUATIONS OF MOTION 



Axial Force Equation 

• • • 2 2 

mu - my r + mz„q = mvr - mwq + mx„ (q + r ) - my pq - mz„pr 

G Vj G G G 

+ ^ i^fx' q^ + X’ + X' rpl 
2 |_ qq^ rr rp 

+ ^2,^fx-'u + X' vr + X' wql 
2 L u vr wq J • 

+ J i^\x' X? + X' + X' w^l 

2 L uu vv ww J 

+ Y l^thrcst ^a^ + b^n + 

+ ^ l^{r\ - l)Fx' u^(Dl^ + D2^ + D3^) + X’ 

2 •^L 'S6q wn 

- (W - B) sin 0 - drag 



+ X' w 
wwq 
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« 






hO|TD hOj-O t^[0 



Lateral Force Equation 



mu - 



mz^p + mx^r = 



2 2 

mwp - mur + my^(r + p ) - mz^gr - mx^qp 

+ ^ a^fY* 'p +Y*'r+Y'| iplpl +Y' pq+Y' qrl 
2 Lp*^ r pIpT'^^' pq*^^ qr^ J 

+ ^ 2,^1 Y* 'v + Y' vq + Y' wp + Y' wr + Y'ur + Y'up 
2[_v vq wp^ wr r p^ 

+ Yj^|^u|r| (D1 - 0.6088D2 + 0.6088D3) 



+ ^ (v^ + w^)^|r| + Y;^ur(n 



rn 

1 



->] 



^ il^jY^u^ + Y^uv + Y_'j_|V(v^ + 



V V 



+ Y' 
vw 



vwu 



2 2 




U + V 


+ 1 


Y' uv + 


Y’ 


. vn 


V 


+ Y' 

5ri 


(D1 



+ Ylu^CDl - 0.6088D2 + 0 
0 



2,2 



V V n 



+ (W - B) sin (() cos 0 



ww 



/ D(x) dx 

a i 



J H(x)v(x)| 



[w(x)]^ + [v(x)]^p dx 



X 



w 



fw 



(x^ - X / 

fw so X 



>/ 



w(x)v^^(t - t(x)) dx 



vs 



.6088D3) 



88 



hjp 



Normal Force Equation 



niw - 



mx^q + my^p = 



2 2 

muq - mvp + mz (p + q ) - nix^rp - my_rq 

Cjr (j 



+ I + z; 



[ Z.’a + 7' P^U 

n/u + V + w 

+ z- 1 

2^2 

V U + V + w J 



+ Z’ I , 

r r 



2 

r u 



./2 2 ^ 2 
V U + V + w 



+ 2 ^ 



z* ’w + z’ 
w vr 



vru 



/ 2 ^ 2 ^ 2 
N/u + V + W 



+ Z’ vp + Z’uq 
vp ^ q 



+ Zj^|^u|q| (-0.7934D2 - 0.7934D3) 



+ Z'l I (v^ + w^)^|q| + Z’ uq(n 
w|q||w|^ qn 






+ 2 ^ 



Z!u^ + Z'uw + Z' I |w(v^ + + Z| lulwl 

* w w w w ' ' 



+ Z' I I 

V V 



2 

V u 



+ Zlu^(-0.7934D2 - 0.7934D3) 



J 2 ^ 2 ^ 2 

WU + V + w 



+ Z' u^6, + Z’ |w 
ob b ww' 






^ (n - 1) 



2 2 2 

Z'uw + Z'i |W(v +w) 
wn w w q 



+ Zl u^(-0.7934D2 - 0.7934D3)! + (w - B) cos i 
on 



>] 



cos 4> 
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h 0|-0 



J C^D(x)w(x) |[v(x) + [w(x)]^p dx 

£ 



£Cv J v(x)v^^(t - Z^) dx 

X 

vs 
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Rolling Moment Equation 



my^w 



mzcv = 



~I p — (I — I ) gr + (r + pq)I — (r^ - q^)I — (pr — q)I 
xx^ zz yy xz ^ ^ yz xy 



+ my^uq - my^vp - mz^wp + mz^ur 
+ T K* 'P + K- ’r + K' qr + K' pq + K' , iplpl 

2 [P P '1 pIpT J 

+ ^ Jl^lK'up + K'ur + K* 'v + K' vq + K' wp + K' wr 
ZLp'^ r V vq^ wp^^ wrj 






Kj.u^ + K'uv + Kj^|^jv(v^ + + Kj_,vw + K{u^(Dl - D2 + D3) 



V V 



vw 



+ K^^u^Cn - 1)^ - Propt + 

(k 4S6^ sin 4((t> - 37.5) + K8S sin 8((}) - 37.5)) 

\ s s s / 

+ [ (y^w - yTjB) cos 0 cos cj) - (z^w - z B) cos 0 sin cj)] 

(j D (j D 

_^fw 

+ ^ l^C z f w(x)v^ (t - T ) dx 
2 J fw x^ 

X 

vs 

+ I + ^l) 

sin (4(<^g - 37.5°)) + sin (8(c(>^ - 37.5°))] 



V = V + X r 
s s 



w = w - X q 
s s^ 
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6 

s 



-1 

tan 




+ w 



f 



u 



I 



I 
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Pitching Moment Equation 



mzcu 



mxgw = 



mz^vr - mz^wq - mx^uq + mx^vp - I^^q - + I^^(p + qr) 

- 



+ 2-1^ 



M- 'q + (m’ I IP^ + M' I ir^'i 

q ^ \ pIpI r|r| / 



J 2 , 2 ^ 2 

N/U + V + W 



+ M^^rp + |^|q| q 



rp 






M-'w + M' +M + M’uq 

2 , 2 ^ 2 ^^J~2—r—2 ^ 

Vu+v+w N/u+vH-W 

+ Mj^j^^ujql (Dl^ - 0.6088D2^ - 0.6088D3^) 



+ Mj^j^u|q| (0.7934D2 + 0.7934D3) + Mj^|^q(v^ + 
+ M' uq(n - 1) 

qn 



+ I 



M'u^ + M'uw + M' I |w(v^ + + mI iuIwI 

* w w w w ' ' 



+ M' I 

V V 



+ Mlu^(0.7934D2 + 0.7934D3) + M’ u"5, + M’ 

0 Ob D w 



V u 



J 2 2 ^ 2 

V U + V + w 



+ M' u^(Dl^ - 0.6088D2^ - 0.6088D3^) 
00 



. .2, 



ww 



w|Jv^ + w^j 



|w 



J!,^(ri - 1 )|m' ijw + M' I I w(v^ + 

' wn w w n 



+ M' u^(Dl^ - 0.6088D2^ - 0.6088D3^) 
66n 



+ M^^U^(0.7934D2 + 0.7934D3) 
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top top 



(x^w - ^ ^ ^ 



C^ /c^xD(x)w(x)|[w(x) ]^ + [v(x)]^p dx 



fw 

iC f xv(x)v^ (t “ T ) dx 
V J fw X 



X 



vs 
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ro|o 



Yawing Moment Equation 



mx^v - my^u = 



mx wp - mx ur - my vr + my wq - I r-(I -I )pq + I (q + rp) 



zz 



yy XX' 



yz 



-I (q^ - p^) - I (rq - p) 
xy ^ ^ xq ^ ^ 



+ ^ n’ 'r + N* 'p + N’ pq + N' qr + N ' i irlrll 
2 L ^ P pq r|r| ‘ j 

+ ^£^ N*’v+N’ wr+N’ wp+N’ vq+N’up+N’ur 
2Lv wr wp^ vq p^ r 

+ Nj^|^^u|r| (0.7934D2^ - 0.7934D3^) 

+ n! |^u|r| (D1 - 0.6088D2 + 0.6088D3) + n! , r(v^ + 
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APPENDIX C 



LINEARIZED EQUATIONS OF MOTION 
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Normal Force 
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+ 2“o“s®10 - O-SOSSDjj, + 0.608803^) + - 1) 

2V55„(i„ - 1)(0.7934D^„ - 0.7934d2J 

- l)(Dj(, - O. 6 O 88 D 3 J + 0 . 6088 D 3 „)] 



+ Av 



2 2 

r V 2v + w 

00 00 

“iny^r + N q + N i i + N u + N i i 

•^G o vq^o |v|r v o vlv| 



/T “""2 

/v + w 
^ o o 



/ 2 ^ 2 
^ o o 



2^2 
u w u + w 

+ N — — — + N u (n - 1) + N I I 

vw _3 vqo o v|v|ri 

2\2 



/ 2 ^ 2 ^ 2 \ 
U + V + w 

\ O O 0/ 



O 2 2 

2v + w 
o o 



/ 2 ^ 2 
Vo o 



8 N 



CRFW 



3v 



+ Aw 



mx^p + iny^q +N r +N p +Ni i 
G^o wr o wp^o v r 



r w 
o o 



/ 2 ^ 2 
. /V + w 
Vo o 



+ N 



V w 
o o 



V V 



I 2 ^ 2 

^/v + w 

Vo o 



+ N 



/ 2 2 \ 

U V u + V 
o o \ O Of 



vw 



3 
2\2 



/ 2 ^ 2 ^ 2 \ 
U + V + w I 
\ o o O' 
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« 



I 

! 

I 



- «'' o”o . ®''CRFW 

+ N I I H r 

vlvjri /-^ 2 

V + w 
o o 



+ Ap[mx„w - (I -I)q+Ir+2Ip+Nq+Nw 
Go yy XX yz o xy o pq o wp o 



+ Aq 



[" yc"o - 



(I -I)p-2Iq-Ir+Np+Nr 
yy XX o xy^o xy o pq o qr o 



+ N V + 



vq o 3q 



^^CRFwl 
9q J 



+ Ar 



-mx-u - my_v + I p - I q +N q +2Ni ilrl+N w 
Go *^G o yz^o xz o qr^o ^1^1 ° wr o 



+ •'r-'o * X|r|6 S tIt 



•|66 



"irU tH ®10 • “■“ 8®20 + 0 -“ 88D3 „) + 

II o 

+ N . 1 ) 

rq o' o 9r J 

+ A<i)[(XgW - XgB) cos 0^ cos <()^] 

+ A0[-(XgW - XgB) sin 0^ sin + (YqW - YgB) cos 0^] 

+ “ 6 “o *' 6 n “3 

A«2[i-5868N|^|jjU^|t^|D 2„ - 0.6088N | ^ | juJ r^| + l.SSASNjjU^Dj^, 

? 2 2l 

- 0.6088N^u, + 1.5868N^^„u„D^^ - 0.6088N^^u^J 



+ A6 



'6 o 



' 66 n o 20 



*«3[-l-5868N|^|jjV.JrJ + 0.6088N| ^ | 

+ O.6O88N5 UJ - 1-5868Nj5^uJd 3„ + O.AOSSNj^u^] 



r - 1.5868N.;^u D_. 
o' 00 o 30 
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'^''["rn"o'o * ''vn“o''o "v|v|n''J''o ”i 

+ "«Sn"o(»-^”4D2 . 0.7934 d2 ) 



+ N. u (D-^ - 0.6088D^^ + 0.6088D_) 
6n o 10 20 30 , 



where: 



3N 



CRFW 



3v 



8N 



CRFW 



8w 



(“ 



ww 



/ D(x) dx 
i 



= K- 



ww 

/ D(x) dx 
Z 



)( 

)/ 



-f 2v^(!;)) 



dx 



V V (x)w (x) 
X H(x) o o 



JvJ(x) +wj(x) 



dx 



8N 



CRFW 



8q 



9N 



CRFW 



8r 






ww 



i 



2 , V (x)w (x) 

^ T^/ N , I f X H(x) o o j 

/ D(x) dx/7 dx 



= f(- 



WW 



/ D(x) dx 



)/ 



X H(x) 



n/v^x) + wj(x) 

(w^(x) + 2v^(x)) 
JvJ(x) + w^(x) 



dx 
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APPENDIX D 



SUMMARY OF NON -DIMEN SIGNALIZING FACTORS AND 
DIMENSIONALIZED HYDRODYNAMIC COEFFICIENTS 



NOTATION [1] 



Symbol Dimensionless Form 



Definition 



I’ 

X 



X 



I* 

y 



ip*' 



I' 

z 



ip*' 



xy 



I’ 

xy 



J2L 






yz 



I’ 

yz 






ip*' 



zx 



I’ 

zx 



zx 



ip*' 






k; = 






1 3 2 



Moment of inertia of submarine about 
X axis 



Moment of inertia of submarine about 
y axis 



Moment of inertia of submarine about 
z axis 



Product of inertia about xy axis 



Product of inertia about yz axes 



Product of inertia about zx axes 



Rolling moment when body angle (a, 3) 
control surface angles are zero 



and 
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NOTATION (Cont.) 



Symbol 


Dimensionless Form 


Definition 


K* 

*n 




j pjiV 


Coefficient used in representing as a 

function of (n-1) 


K 

P 


K' = 
P 


K 

. P 


First order coefficient used in represent 
ing K as a function of p 


j PZ^U 


K* 

P 


K* ’ = 

P 


K* 

P 


Coefficient used in representing K as a 
function of p 


""pIpI 


K'l 1 

PlPl 


. !eM 

Y pz^ 


Second order coefficient used in repre- 
senting K as a function of p 


K 

pq 


K' 

pq 


K 

pq 


Coefficient used in representing K as a 
function of the product pq 


K 

qr 


K' 

qr 


■'pr 

ip^= 


Coefficient used in representing K as a 
function of the product qr 


K 

r 


K' = 
r 


K 

r 


First order coefficient used in represent- 
ing K as a function of r 


j pz^u 


K- 

r 


K- ' = 

r 


K* 


Coefficient used in representing K as a 
function of f 


K 

V 


K' = 

V 


K 

V 


First order coefficient used in represent- 
ing K as a function of v 


j 



i) 

I 

ii' 



r 



4 



NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 



K- 



V 




Coefficient used in representing K as a 
function of v 



K 



V V 



K’ 



V V 



K 



V V 




Second order coefficient used in repre- 
senting K as a function of v 




K’ 



vq 




Coefficient used in representing K as a 
function of the product vq 



K 



vw 



K’ 



vw 




Coefficient used in representing K as a 
function of the product vw 



K 

K K' = 7 

wp wp 1 



Coefficient used in representing K as a 
function of the product wp 



K 



wr 



K' 



wr 



K 



wr 




Coefficient used in representing K as a 
function of the product wr 



^6 



^6 



K 



6r 



1 3 2 

j pil V 



First order coefficient used in represent- 
ing K as a function of 6^ 



Z Z^ = 1 



Overall length of submarine 



m 



m 



m 




Mass of submarine, including water in 
free-flooding spaces 



M. 









3 2 
pZ V 



Pitching moment when body angles (a, 6) 
and control surface angles are zero 
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NOTATION (Cont.) 



Syn±>ol 


Dimensionless Form 


Definition 


M 1 1 

pIpI 


M'l 1 

pIpI 


= 


Second order coefficient used in repre- 
senting M as a function of p. First order 
coefficient is zero. 


M 

q 


M’ = 

q 


M 

q 


First order coefficient used in represent- 
ing M as a function of q 




M- 

q 


M- ' = 

q 


M- 

a_ 


Coefficient used in representing M as a 
function of q 


“q|q| 


""qkl 


Y Pi^ 


Second order coefficient used in repre- 
senting M as a function of q 




""Iqk 


- ""kk 

Y 


Coefficient used in representing as a 

function q 


M 

rp 


M' 

rp 


M 

Y p^^ 


Coefficient used in representing M as a 
function of the product rp 


M 1 1 
r|r| 


M' 1 1 
r|r| 


- ""rlrl 

■ 


Second order coefficient used in repre- 
senting M as a function of r. First order 
coefficient is zero 


M 

vr 


M' 

vr 


M 

vr 


Coefficient used in representing M as a 
function of the product vr 


M 1 1 
v|v| 


M' 1 1 
v|v| 


- "v|v| 


Second order coefficient used in repre- 
senting M as a function of v 



I 



s 

t 

t 

\ 

< 



NOTATION (Cont.) 



Symbol 



M 

w 



Dimensionless Form 
M 



M’ = 
w 



w 






Definition 



First order coefficient used in represent- 
ing M as a function of w 



M 



wp 



M* 

wri 



M 



wn 



j P£^U 



First order coefficient used in represent- 
ing as a function of (p - 1) 



M* 



w 



w 



M- 



w 






Coefficient used in representing M as a 
function of w 



Ml 



w 



“(»i 



Ml 



w 






First order coefficient used in represent- 
ing M as a function of w; equal to zero 
for symmetrical function 



Mi 



|w|q 



Mf 



|w|q 



Ml 



|w|q 






Coefficient used in representing M as a 
function of w ^ 



M 



w w 



M’ I I 

w w 



M , , 

w| w| 

I pa^ 



Second order coefficient used in repre- 
senting M as a function of w 



M 



w w n 



M’ I I 

w w n 



M I ( 
w| w| p 

2 P ^ 



First order coefficient used in represent- 
ing Ml I as a function of (p - 1) 

® w w 



M 



6b 



M 



^6b 



5b 



1 3 2 

2 U 



First order coefficient used in represent- 
ing M as a function of 6^ 



M, 



= 

^ ^ p.V 



First order coefficient used in represent- 
ing M as a function of 6 
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NOTATION (Cont.) 



Symbol Dimensionless Form 



M 



6n 



MI 



6n 



1 3 2 

j Pi^ V 



Definition 



First order coefficient used in represent- 
ing M^ as a function of (n - 1) 



N. 



, 

= 1 TT 

J pH V 



Yawing moment when body angles (a, 6) and 
control surface angles are zero 



N 



N 



N» = c 

p 1 4 



First order coefficient used in represent- 
ing N as a function of p 



N- 



N- 

N- ' = P 



Coefficient used in representing N as a 
function of p 



N 



pq 



N 



N’ 

pq 



pq 






Coefficient used in representing N as a 
function of the product pq 



N 






N 



N’ 

qr 






1 ,5 

j 



Coefficient used in representing N as a 
of the product qr 



N 



N 



N’ = ^ ^ 

r 14 



First order coefficient used in represent- 
ing N as a function of r 



N 



rn 



N’ 

rn 



rn 



1 

J pZ V 



First order coefficient used in represent- 
ing N^ as a function of (n “ 1) 



N- 



N- 

N.' = 

r 

2 P" 



Coefficient used in representing N as a 
function of f 
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J 



! 



! 



NOTATION (Cont,) 



Symbol Dimensionless Form 



Definition 



u 

u 


N’ 1 1 
r 1 r 1 


= -L[r| 


Second order coefficient used in repre- 
senting N as a function of r 


1 r 1 6r 


1 r 1 6r 


- ""khr 


Coefficient 
function of 


used in representing N^^ as a 


N 

V 


N’ = 

V 


N 

V 


First order coefficient used in represent- 
ing N as a function of v 


Y 


N 

vri 


N» 

vn 


N 

vn 

Y 


First order 
ing N as a 


coefficient used in represent- 
function of (n - 1) 


N- 

V 


N- ’ = 

V 


N* 

V 

1 


Coefficient 
function of 


used in representing N as a 

• 

V 


N 

vq 


N» 

vq 


N 

vq 


Coefficient 
function of 


used in representing N as a 
the product vq 


N| 1 
1 v| r 


Nf 1 
|vlr 


. !kk 


Coefficient 
function of 


used in representing N^ as a 

V 



N , 

V I V 



N», 
v| V 



N 



V V 




Second order coefficient used in repre- 
senting N as a function of v 



N 1 I 
V 1 v| n 



N’l , 
v| v|n 



N I I 
V V h 




First order coefficient used in represent- 
ing N 1 I as a function of (n - 1) 

® v| v| 
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NOTATION (Cont.) 



Symbol Dimensionless Form Definition 



N 

vw 



N 

vw 


N' 

vw 


vw 

1 .3 

2 P" 


Coefficient used in representing N as a 
function of the product vw 


N 

wp 


N’ 

wp 


N 

wp 


Coefficient used in representing N as a 
function of the product wp 


N 

wr 


N’ 

wr 


N 

wr 


Coefficient used in representing N as a 
function of the product wr 




N! 

5r 


“jr 

i P*V 


First order coefficient used in represent 
ing N as a function of 6^ 


^6rn 


N ' = 


1 3 2 

2 P«^ U 


First order coefficient used in represent 
ing as a function of (n - 1) 


P 


P’ = 


pi 

U 


Angular velocity component about x axis 
relative to fluid (roll) 


• 

P 


• 

p’ = 




Angular acceleration component about x 
axis relative to fluid 


q 


q’ = 


u 


Angular velocity component about y axis 
relative to fluid (pitch) 


• 

q 


• 

q' =- 


qS? 

u' 


Angular acceleration component about 
y axis relative to fluid 


r 


r ’ = 


rit 

U 


Angular velocity component about z axis 
relative to fluid (yaw) 



123 




I 



i 

I 

> 

i 
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NOTATION (Cont.) 



Form 



Definition 



Angular acceleration component about 
z axis relative to fluid 



Linear velocity of origin of body axes 
relative to fluid 

Component of U in direction of the x axis 



Time rate of change of u in direction of 
the X axis 



Command speed: steady value of ahead 

speed component u for a given propeller 
rpm when body angles (a, 3) and control 
surface angles are zero. Sign changes 
with propeller reversal. 



Component of U in direction of the y axis 

Time rate of change of v in direction of 
the y axis 

Component of U in direction of the z axis 

Time rate of change of w in direction of 
the Z axis 

Longitudinal body axis; also the coordi- 
nate of a point relative to the origin of 
body axes 



The X coordinate of CB 






I 



I 



NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 




The X coordinate of CG 



X 

- - 2 . 



qq 



X* 

qq 



qq 



rp 



X’ 

rp 



X 

1 

2 



rr 



X' 

rr 



rr 



X- 



x-» 

u 






UU 



X' 

UU 



X 

UU 

2 



A coordinate of the displacement of CG 
relative to the origin of a set of fixed 
axes 



Second order coefficient used in repre- 
senting X as a function of q. First order 
coefficient is zero 



Coefficient used in representing X as a 
function of the product rp 



Second order coefficient used in repre- 
senting X as a function of r. First order 
coefficient is zero. 



Coefficient used in representing X as a 
function of u 



Second order coefficient used in repre- 
senting X as a function of u in the non- 
propelled case. First order coefficient 
is zero . 



X 



vr 



X’ 



vr 



X 

vr 




Coefficient used in representing X as a 
function of the product vr 



X 



w 



X’ 



vv 



X 

w 




Second order coefficient used in repre- 
senting X as a function of v. First order 
coefficient is zero. 
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NOTATION (Cont,) 



Symbol Dimensionless Form 



Definition 



X 



vvri 



X’ 

wn 



wn 




First order coefficient used in represent- 
ing X^ as a function of (n - 1) 



X 

X X* = 

wq wq l„o3 

2 



Coefficient used in representing X as a 
function of the product wq 



X 



ww 



X’ 

ww 



X 

ww 




Second order coefficient used in repre- 
senting X as a function of w. First order 
coefficient is zero 



X 



wwn 



X’ 

wwn 



wwn 




First order coefficient used in represent- 
ing X as a function of (n - 1) 

® ww 



^6b6b 



Y » 

6b5b 



X 



6b5b 



1 2 2 



Second order coefficient used in repre- 
senting X as a function of 5^. First 
order coefficient is zero. 



6r5r 



XI . 

6r6r 



X 



6r6r 



1 .2^2 
“2 U 



Second order coefficient used in repre- 
senting X as a function of 6^. First 
order coefficient is zero. 



Xr X K z 
ororn ororn 



5r5rn 
1 2 2 



First order coefficient used in represent- 
ing X. . as a function of (n - 1) 

^ 5r5r 



Xx K z 

osos 5s6s 



X 



6s5s 



1 2 2 
Y pil V 



Second order coefficient used in repre- 
senting X as a function of 6^. First 
order coefficient is zero. 



X 



6s6sn 



ososn 



5s5sn 

1 2 2 
j Pi^ U 



First order coefficient used in represent- 
ing X^ ^ as a function of (n “ O 
^ 6s6s 
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NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 



y 



y’ 



= Z 



i 



Lateral body axis; also the coordinate of 
a point relative to the origin of body 
axes 



y 



B 




The y coordinate of CB 







The y coordinate of CG 



y 



o 



Y 



* 



Y 

P 



Y- 

P 



y^ = ^ coordinate of the displacement of CG 

relative to the origin of a set of fixed 
axes 



Y ’ 



Y’ 

P 



Y- ’ 
P 





Y 


1 


.2„2 


2 


pS, U 




Y 




P 


1 


3 


2 


p£ U 




Y* 




P 



Lateral force when body angles (a, B) and 
control surface angles are zero 



First order coefficient used in represent- 
ing Y as a function of p 



Coefficient used in representing Y as a 
function of p 



Y 





Second order coefficient used in repre- 
senting Y as a function of p 



Y 

pq 



Y’ 

pq 




Coefficient used in representing Y as a 
function of the product pq 
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f 



NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 



Y 

qr 




Coefficient used in representing Y as a 
function of the product qr 



Y 



r 



Y’ 



r 



Y 

r 






First order coefficient used in represent- 
ing Y as a function of r 



Y 

rri 




Y 

rn 






First order coefficient used in represent- 
ing Y^ as a function of (n - 1) 



Y- 

r 




Coefficient used in representing Y as a 
function of f 



r |6r 






U 



Coefficient used in representing Y^^ as a 
function of r 



Y 

V 



Y’ 



V 



Y 

V 



Y pii u 



First order coefficient used in represent- 
ing Y as a function of v 



Y 

vn 




Y 

vn 



1 

Y p£ u 



First order coefficient used in represent- 
ing Y^ as a function of (n - 1) 



Y- 

V 




Coefficient used in representing Y as a 
function of v 



Y 

vq 



Y’ 

vq 




Coefficient used in representing Y as a 
function of the product vq 
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NOTATION (ContO 



Symbol 



V r 



Dimensionless Form 
Y 






v|r I 

ip*" 



Definition 



Coefficient used in representing Y as a 
function of r ^ 



V V 



nivi 



^i''i 



Second order coefficient used in repre- 
senting Y as a function of v 



V V n 



Y’ , 1 
V V n 



v| vln 

ip*" 



First order coefficient used in represent- 
ing ^vlvl ^ function of (n - 1) 



vw 



Y’ 

vw 



vw 



Y 09 } 



Coefficient used in representing Y as a 
function of the product vw 



wp 



Y’ 

wp 



wp 



ip*" 



Coefficient used in representing Y as a 
function of the product wp 



wr 



Y’ 

wr 



wr 



1 n3 

2 



Coefficient used in representing Y as a 
function of the product wr 



■6r 



or 



■6r 



1 2 2 
y U 



First order coefficient used in represent- 
ing Y as a fiinction of 6r 



“6rri 



orn 



6rri 



1 . 2,,2 
Y U 



First order coefficient used in presenting 

Y^ as a function of (n - 1) 

6r 



Normal body axis; also the coordinate of a 
point relative to the origin of body axes 
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I 
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NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 



z 



B 




The z coordinate of CB 



z 



G 




The z coordinate of CG 



z 

o 




A coordinate of the displacement of CG 
relative to the origin of a set of fixed 
axes 



Z 



* 



z: 

* 



1 2 2 

J pi 



Normal force when body angles (a, g) and 
control surface angles are zero 





Second order coefficient used in repre- 
senting Z as a function of p. First order 
coefficient is zero 



Z 

q 



z’ 

q 



z 

q 






First order coefficient used in represent- 
ing Z as a fianction of q 



Z- 

q 




Coefficient used in representing Z as a 
function of q 



lql6s 



|ql6s 



lql<5s 



Coefficient used in representing Z^^ as a 
function of q 



Z 

rp 





Coefficient used in representing Z as a 
function of the product rp 
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NOTATION (Cont.) 



Symbol 


Dimensionless Form 


Definition 


z 1 1 

r|rl 


^r|rl 


= fikL 


Second order coefficient used in repre- 
senting Z as a function of r. First order 
coefficient is zero 


z 

w 


z’ = 

w 


Z 

w 


First order coefficient used in represent- 
ing Z as a function of w 


1 

Y pil u 


Z 

qn 


z* 

qn 


z 

an_ 

1 PJI^U 


First order coefficient used in represent- 
ing Z^ as a function of (q - 1) 


z 

W71 


z’ 

wq 


z 

wq 

■J 


First order coefficient used in represent- 
ing Z as a function of (q - 1) 
w 


z- 

w 


Z* ’ = 

w 


Z- 

w 


Coefficient used in representing Z as a 
function of w 


^|w| 


^|w| 


= 

y pJl^U 


First order coefficient used in represent- 
ing Z as a function of w; equal to zero 
for symmetrical function 


^w|q| 


^w|q| 


= ^w|q| 


Coefficient used in representing Z^ as a 
function of q 


^w| w| 


^w|w| 


^w| w| 

1 p2 
2 


Second order coefficient used in repre- 
senting Z as a function of w 


1 w 1 n 


Z' 1 1 

w| w| q 


Z 1 1 

w| w| q 

1 p2 

IP' 


First order coefficient used in represent- 
ing Z 1 I as a function of (q - 1) 
w w 



NOTATION (Cont.) 



Symbol Dimensionless Form 



Definition 



Z 

ww 



Z 



6b 



Z 



6s 



Z’ 

ww 




Second order coefficient used in repre- 
senting Z as a function of w; equal to 
zero for symmetrical function 



Z ’ 

^6b 



"6b 



1 2 2 

Y pii V 



First order coefficient used in represent- 
ing Z as a function of 6, 

D 




"6s 



1 0 2,,2 
Y p£ u 



First order coefficient used in represent- 
ing Z as a function of 6^ 



^6sn 



6sn 



"6sn 



1 2 2 
2 U 



First order coefficient used in represent- 
ing Z. as a function of (n - 1) 
os 



6 



b 



Deflection of bowplane or sailplane 



6 Deflection of rudder 

r 



6 Deflection of sternplane 

s 



n 



The ratio 



u 




0 



Angle of pitch 






Angle of yaw 



<P 



Angle of roll 



a . 

1 






Sets of constants used in the representa- 
tion of propeller thrust in the axial 
equation 
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Axial Coefficients 



XQQ 


XRR 


XRP 


XUDOT 


-138568. 


-11547. 


54588. 


-1819.584 


XVR 


XWQ 


XUU 


XW 


63161. 


-58728. 


0 


154.872 


xww 


XDD 


XWN 


XWWN 


120.204 


-188.02 


0 


0 


XDDN 


XDBDB 






0 


-91.044 







Lateral Coefficients 



YRDOT 


YPDOT 


YPABP 


YQR 


-39890. 


-72433. 


0 


0 


YVDOT 


YVQ 


YWP 


YWR 


-79122. 


0 


58728. 


0 


YR 


YP 


YABRD 


YVABR 


41705. 


-4549. 


0 


0 


YSTAR 


YV 


YVABV 


YVW 


0 


-1043. 


-2054.2 


0 


YD 


YRN 


YVN 


YVABVN 


185.036 


0 


0 


0 


YDN 


YPQ 






43.371 


138568. 







133 



Normal Coefficients 



ZQDOT 


ZPABP 


ZRABR 


ZRP 


-170000. 


0 


0 


0 


ZWDOT 


ZVR 


ZVP 


ZQ 


-69459. 


0 


-33067. 


-42148. 


ZABQDD 


ZWABQ 


ZSTAR 


ZW 


0 


0 


0 


-1044.3 


ZWABW 


ZABW 


ZVABV 


ZD 


0 


-8.97 


0 


248.054 


ZQN 


ZWN 


ZWABWN 


ZDN 


0 


0 


0 


31.34 


ZDB 


ZABQD 


zww 




-175.90 


0 


15.99 






Rolling Coefficients 




KPDOT 


KRDOT 


KQR 


KPQ 


1133000. 


-906993. 


-2796560, 


0 


KPABP 


KP 


KR 


KVDOT 


0 


-153265. 


-398900. 


-15746. 


KVQ 


KWP 


KWR 


KSTAR 


0 


53283. 


0 


0 


KV 


KVABV 


KVW 


KD 


-3423. 


0 


0 


1568. 


KSTARN 


K4S 


K8S 




0 


1281.4 


0 
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Pitching Coefficients 



MQDOT 


MPABP 


MRABR 


MRP 


-124300000. 


0 


0 


29666000. 


MQABQ 


MWDOT 


MVR 


MWW 


0 


-170600. 


0 


-1595.6 


MQ 


MABQDD 


MABQD 


MABWQ 


-4066000. 


0 


0 


0 


MSTAR 


MW 


MWABW 


MABW 


0 


221.62 


-35873. 


-239.72 


MVABV 


MDD 


MD 


MQN 


0 


0 


-23438. 


0 


MWN 


MWABWN 


MDDN 


• MDN 


0 


0 


0 


-3058. 



MDB 

4222. 



Yawing Coefficients 



NRDOT 


NPDOT 


-128800000. 


-434600 


NRABR 


NVDOT 


0 


-39800. 


NVQ 


NP 


0 


-62986, 


NABRD 


NABVR 


0 


0 


NVABV 


NVW 


49980. 


0 


NRN 


NVN 


0 


0 



NDN 

-4374 

CDCRSS 

.56889 



NPQ 


NQR 


-97313000. 


0 


NWR 


NWP 


0 


-138558, 


NR 


NABRDD 


-3766000. 


0 


NSTAR 


NV 


0 


-12661. 


NDD 


ND 


0 


-9798. 


NVABVN 


NDDN 


0 


0 
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APPENDIX E 



DESCRIPTION OF PROPULSION MODELS AND NOMINAL POINT DETERMINATION 



There are two propulsion models an RPS and an eta propulsion 

model. The RPS model is a first order differential equation in terms of 

RPS (revolutions per second) and is considered the more accurate model 

of the two. The eta propulsion model is a simplified version of the RPS 

model in that it is a first order differential representation where eta 
u 

(n) is defined as — . U is the actual speed of the submersible and u^ 
is the commanded forward velocity. The eta model is important because 
it was the model that was linearized and included in the linear equa- 
tions of motion. 

In order to find a nominal point, the final values of the states 
(u, V, w, p, q, r, phi and theta) from the RPS integration are used as 
an initial point with the eta model integration. The nonlinear equa- 
tions are again integrated, however, using the eta propulsion model. A 
search routine is used to search for a nominal point (a local equilib- 
rium point where the derivatives of the state variables are zero) . To 
get some idea how good the nominal point is, the accelerations asso- 
ciated with each of the states -are calculated. As a rule of thumb, 
accelerations should be less than 10 . The equilibrium point for this 

model is listed below. 
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XDOT 


(1) 


= 


-0.134572E-10 


XDOT 


(2) 


= 


-0.132110E-07 


XDOT 


(3) 


= 


-0.130859E-08 


XDOT 


(4) 


= 


-0.297533E-08 


XDOT 


(5) 


= 


0.137816E-10 


XDOT 


(6) 


= 


0.999896E-10 


XDOT 


(7) 


= 


-0.568434E-13 


XDOT 


(8) 


= 


0.145519E-10 




u 


= 


0.351355E+02 




V 


= 


0.317284E+00 




w 


= 


-0.6 7850 8E-02 




p 


= 


0.368353E-06 




q 


= 


-0.996159E-04 




r 


= 


-0.368482E-02 


phi 


= 


0.154856E+01 


theta 


= 


0.572548E-02 


psi 


= 


-0.211 14 lE+02 




z 


= 


-0. 170338E-H)0 




DB 


= 


0 




D1 


= 


2.0 




D2 


= 


0.5 




D3 


= 


-0.4744 
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APPENDIX F 



0B8ERVABILITY(T*C) 

-2.1290E-05 7.15B3E-05 
O.OOOOE+00 O.OOOOE+00 

2.6489E-03 -2.0743E-03 
O.OOOOEtOO O.OOOOE+00 

-7.B555E-01 1.6809E-03 
O.OOOOE+00 O.OOOOE+00 



C0MTR0LLABILITY<T"‘«B1 

-1.5311E-02 2.3498E-01 
-8.0762E-11 7.1922E-11 

4.2570E-01 -6.0879E-02 
1.1603E-10 -9.0817E-U 

-3.26B2E+00 2.0472E+00 
-3.6227E-09 2.9746E-09 

1.0639E+00 -7.2540E-02 
1.1672E-10 1.1699E-10 

3.0824E-ai -1.0009E-01 
9.6169E-01 -8.6784E-01 

3.0824E-01 -1.0009E-01 
•9.6169E-01 8.&784E-01 

-1.6724E-01 2.3913E+00 
1.7908E-10 -1.7178E-10 



MODAL ANALYSIS MATRICES 



4.8423E-03 5.2B30E-03 
O.OOOOE+00 O.OOOOE+00 

-1.9729E-02 2.4111E-03 
O.OOOOE^O O.OOOOE+00 

1.3494E-02 3.34SBE-02 
O.OOOOE^OO O.OOOOE+00 



3.2220E-03 3.2220E-03 
B.9118E-03 -B.9118E-03 

-5.BS58E-01 -S.6558E-01 
7.0275E-02 -7.0275E-02 

1.4757E-03 1.4757E-03 
3.1573E-04 -3.1573E-04 



1.3762E-03 -2.7719E-02 
O.OOOOE+00 O.OOOOE^OO 

-6.4347E-03 2.6790E-01 
O.OOOOE+OO O.OOOOEm 

6.7821E-02 4.3741E-03 
O.OOOOEtOO O.OOOOE+00 



2.1431E-01 

-7.5965E-11 

7.5270E-02 

1.17&2E-10 

-1.7453E+00 

-3.3771E-09 

1.127iE+00 

-7.2679E-11 

3.1597E-02 

8.3702E-01 

3.1597E-02 

-8.3702E-01 

2.24BSE*00 

1.6932E-10 



1.5552E*00 -i.7186E-01 4.8A29E-01 
-4.403iE-09 3.4853E-09 -3.2005E-09 
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1 

i 



f 

1 

''f 

I'; 

I 

I 



DIAGONAL A HATRII«T'Wt) 

-2.0W3E-02 -4.8182E-U 7.3209E-10 -4.3945E-10 
8.0822E-1S 3.1627E-1S 9.9136E-13 2.S088E-14 

-9.8699E-09 -4.3774E-02 -3.8043E-09 9.8277E-10 
-3.7401E-14 -1.4489E-14 -2.0391E-12 -5.7494E-14 

-1.4485E-09 1.4366E-10 -1.8040E-01 S.4190E-10 
S.6934E-13 2.2166E-13 4.8451E-11 1.5S7SE-12 

-5.2042E-09 -3.5774E-10 4.4974E-09 -1.8432E-01 
-5.9672E-14 -6.8722E-14 2.9312E-13 9.8037E-14 

-1.3921E-09 -3.2740E-10 -1.0737E-09 1.4470E-11 
-7.2923E-10 -3.SS68E-10 1.7830E-10 1.6938E-10 

-1.3999E-09 -3.3429E-I0 -1.0920E-09 3.4241E-11 
7.7433E-10 3.4244E-10 -1.9444E-10 -1.4578E-10 

-1.3947E-09 -4.2439E-10 1.4827E-09 -2.5007E-10 
4.4808E-14 -1.5415E-14 -2.1257E-12 4.2141E-15 

-8.4933E-09 2.2854E-10 -8.58S4E-09 9.8474E-10 
-3.8494E-13 3.7923E-13 2.4495E-11 •1.1123E-12 



EIGENVALUES 

-2.0293E-02 -4.3774E-02 -1.8040E-01 -1.8432E-01 
O.OOOOE^OO O.OOOOE+00 O.OOOOE^ O.OOOOE+00 



-2.I759E-10 -2.9193E-10 -7.0223E-09 -4.9357E-09 
3.I149E-09 -3.1450E-09 9.1370E-13 -4.9420E-12 

-4.5444E-08 -4.5325E-08 -3.1495E-09 1.1828E-07 
4.3310E-09 -4.3132E-<I9 -4.2574E-12 1.2084E-11 

1.3579E-08 1.0251E-08 -4.1731E-09 -1.8558E-09 
-3.4447E-09 4.3284E-09 7.3413E-11 -1.2183E-10 

1.41S2E-08 1.4424E-08 1.054SE-08 -5.S337E-08 
7.0181E-08 -7.0459E-08 -1.3197E-11 -3.9280E-11 

•S.0851E-01 3.4939E-09 8.8439E-10 •1.1038EH)9 
5.2247E-01 -1.7384E-08 1.8838E-09 1.0024E-08 

-8.2S12E-09 -S.08S1E-01 S.8881E-10 3.0042E-09 
2.4480E-08 -5.2247E-01 -1.9305E-09 -4.7210E-09 

-1.1248E-08 -1.1079E-08 -4.2278E-01 2.4002E-08 
-1.8744E-09 1.9342E-09 -4.4594E-13 2.1903E-11 

-8.2833E-09 -1.3044E-08 4.4007E-09 -4.S837E-01 
-4.8472E-10 1.1121E-09 1.1024E-U -2.8917E-10 



-3.0851EH)1 -5.08S1E-01 -4.227BE-01 -4.5837E-01 
5.2247E-01 -5.2247E-01 O.OOOOE+00 O.OOOOE+OO 
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EIGENVECTOR (T) 



2.1339E-02 

O.OMOE+M 

1.&&40E-02 

O.OOOOE^O 

4.1799E-01 

O.OOOOE+CO 

-2.9495E-03 

O.OOOOE+00 

1.5937E-02 

O.OOOOE+00 

-2.1290E-05 

O.OOOOE^OO 

2.&489E-03 

O.OOOOE+00 

-7.8555E-01 

O.OOOOE^O 



-9.9996E-01 

O.OOOOE+00 

-8.3246E-03 

O.OOOOE+00 

-2.7814E-03 

O.OOOOEiOO 

9.6991E-05 

O.OOOOE+00 

-4.4024E-05 

O.OOOOEiOO 

7.1583E-0S 

O.OOOOE^O 

-2.0743E-03 

O.OOOOE+00 

1.4809E-03 

O.OOOOE+00 



-5.3213E-02 

O.OOOOEiOO 

-9,3374EH)1 

O.OOOOE+00 

-3.5312EH)1 

O.OOOOE^O 

3.4084E-03 

O.OOOOEfOO 

-2.2315E-03 

O.OOOOEiOO 

4.8423E-03 

O.OOOOEiOO 

-1.9729E-02 

O.OOOOE+00 

1.3494E-02 

O.OOOOE^OO 



9.0521E-03 

O.OOOOE+00 

-4.4859E-03 

O.OOOOEiOO 

-9.9935E-01 

O.OOOOEiOO 

-3.2107E-04 

O.OOOOE+00 

-4.1744E-03 

0.0000E400 

5.2430E-05 

O.OOOOE+00 

2.4111E-03 

O.OOOOE+00 

3.3454E-02 

O.OOOOEiOO 



2.8894E-03 

-1.5787E-03 

4.4394E-01 

1.8740E-01 

1.4231E-01 

7.5758E-03 

2.5089EH)1 

-3.3123E-01 

1.2570E-03 

S.3841E-04 

3.2220E-03 

4.9118E>03 

•5.4558E-01 

7.0275E-02 

1.4757E-03 

3.1573E-04 



2.8894E-03 

1.5787E-03 

4.4394E-01 

-1.8740E-01 

1.4231E-01 

-7.5758E-03 

2.5089E-01 

3.3123E-01 

1.2570E-03 

-5.3841E-04 

3.2220E-03 

-4.9118E-03 

-5.4558E-01 

-7.0275E-02 

1.4757E-03 

-3.1573E-04 



4.3930E-02 

O.OOOOE+00 

3.9225E-02 

O.OOOOE+00 

9.9S03E-01 

O.OOOOE+00 

4.2573E-03 

O.OOOOE+00 

-4.2192E-02 

O.OOOOE+00 

1.3742EH)3 

O.OOOOEfOO 

-4.4347E-03 

O.OOOOE+00 

4.7821E-02 

O.OOOOE+00 



2.2700E-02 

O.OOOOEiOO 

-9.4447E-01 

O.OOOOEtOO 

4.2817E-04 

0.0000E400 

-1.7434E-01 

O.OOOOEm 

-4.4181E-03 

O.OOOOEiOO 

-2.7719E-02 

O.OOOOE+00 

2.4790E-01 

O.OOOOE^O 

4.3741E-03 

O.OOOOEtOO 
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APPENDIX G 



L MATRIX 



MATRIX 


(14, 3) 

1 


2 


3 


1) 


0.0 


0.0 


0.0 


2) 


0.0 


0.0 


0.0 


3) 


0.0 


0.0 


0.0 




0.0 


0.0 


0.0 


5) 


0.0 


0.0 


0.0 


6) 


l.OOOOOOOOOOOOE+00 


0.0 


0.0 


7) 


0.0 


l.OOOOOOOOOOOOE+00 


0.0 


8) 


0.0 


0.0 


l.OOOOOOOOOOOOE^OO 


9) 


0.0 


0.0 


0.0 


10) 


0.0 


0.0 


0.0 


11) 


0.0 


0.0 


0.0 


12) 


-1 . 835420480000E+01 


-1.403831751000E+00 


-2 . 723467750000E-02 


13) 


-1 . 192239258000E+01 


-1.310703360000E4^00 


-7.691884516000E-02 


1<+) 


1.174827457000E+01 


1.333686037000E4^00 


-4 . 335761088000E-02 
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APPENDIX H 



A, B, C, G, H, (A-BG-HC) Matrices 
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^3703E“01 -.53404E-02 -.19627E-02 0.30208E-02 0.58597E+00 0.57035E400 O.OOOOOE^OO 0.95681E-02 O.OOOOOE^OO O.OOOOOE^OO 
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